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REGIONAL COASTAL PROCESSES NUMERICAL MODELING SYSTEM 


RCPWAVE--A LINEAR WAVE PROPAGATION MODEL 
FOR ENGINEERING USE 


PART I: INTRODUCTION 


1. The "Regional Coastal Processes Numerical Modeling System" research 
work unit is part of the US Army Corps of Engineers (Corps) Shore Protection 
and Restoration research program. Its goal is the development of a modeling 
system that can predict coastal processes on a regional scale so that coastal 
changes resulting from natural forces and man-made structures and modifica- 
tions can be determined on a regional basis. All important physical processes 
that determine coastal changes over a region will be considered in the sys- 
tem's development. Long-wave (e.g. tide and storm surge) hydrodynamics, 
short-wave (e.g. swell and wind sea) propagation, and associated currents and 
setup/setdown will be addressed. These wave phenomena will be used as forcing 
functions to drive models which calculate alongshore and on- offshore sediment 
transport. Sediment sources and sinks such as sediment discharge from rivers, 
dredging gains or losses, and dune erosion will be considered in the sediment 
models. The modeling system will provide a tool for predicting coastal ero- 
sion and deposition, paths of sediment movement, and ultimate fates of coastal 
sediments. 

2. All component models comprising the system must be capable of per- 
forming regional scale simulations in a cost-effective manner. Regional scale 
implies an area of interest with horizontal length scales of up to tens of 
miles and simulation times ranging from days to years. These requirements im- 
pose severe restrictions on candidate models. Individual models must also be 
compatible so that they can interface with one another in an efficient manner. 
Some examples of compatibility are (a) solutions from different models com- 
puted at identical spatial locations and (b) input/output information common 
to different models retrieved/stored using identical formats. 

3. A philosophy for development of the modeling system has been adopted 
which provides direction to the work unit research. The aim is to develop 
(and/or obtain) and link a series of models so that all important physical 


processes are simulated to a predetermined level of sophistication. This 


milestone will occur approximately 2 years before the scheduled conclusion of 
the work unit. Component models comprising the initial system are determined 
by: (a) available models within the Coastal Engineering Research Center 
(CERC) addressing the important physical processes, (b) other (outside CERC) 
state-of-the-art methods for modeling these processes on a regional scale, and 
(c) the need for having a field oriented product in the near future. The 
framework designed for linking these models into a usable system will allow 
for integration of more sophisticated models when they become available. Re- 
search addressing perceived weaknesses in the system models will be conducted 
concurrently with work done to bring the system to its initial level of 
sophistication. 

4, Any technology created as an interim product, which can immediately 
aid Corps field engineers and can be easily transferred to them, will be made 
available. This report documents such a product. The Regional Coastal Pro- 
cesses Wave (RCPWAVE) Propagation Model can be used to predict linear, plane 
wave propagation over a "regional" area of arbitrary bathymetry. This model 
currently forms the initial level of sophistication for the short-wave model- 


ing component of the regional system. 


PART II: REGIONAL COASTAL PROCESSES WAVE PROPAGATION MODEL (RCPWAVE) 


Background Information 


5. Linear wave theory was chosen as the initial level of sophistication 
for the short-wave modeling component because, historically, it has been shown 
to yield fairly accurate first order solutions to wave propagation problems. 
Considering both accuracy and cost, it is currently the most feasible way to 
model waves on a regional scale. Modeling short-wave processes using either a 
fully two-dimensional nonlinear wave theory or a two-dimensional spectral re- 
presentation of irregular waves is presently impractical for the types of ap- 
plications anticipated for this modeling system. 

6. Much of the early work addressing the problem of linear, monochro- 
matic wave propagation was based on wave ray methods and the manual construc- 
tion of refraction diagrams (see Johnson, O'Brien, and Issacs (1948); Dunham 
(1951); and Pierson, Neumann, and James (1952) for examples). During the 
1960's and early 1970's the refraction problem was solved in a more efficient 
way through the use of the computer (for examples see Harrison and Wilson 
(1964), Dobson (1967), Noda et al. (1974), and Rabe (1975). Refraction theory 
fails in regions of complex bathymetry where waves are strongly convergent or 
divergent. Crossing wave rays results in the computation of erroneously large 
wave height estimates. Strongly divergent wave fields manifest themselves in 
regions of unusually small wave heights. Laboratory and prototype observa- 
tions show that refraction theory is inadequate under these conditions (Whalin 
1971 and 1972). 

7. Inclusion of diffractive effects into the equations governing wave 
propagation allows wave energy to be diffused from regions of convergence to 
regions of divergence. Berkhoff (1972 and 1976) derived an elliptic equation 
approximating the complete wave transformation process for linear waves over 
an arbitrary bathymetry constrained only to have mild bottom slopes (hence the 
designation "mild slope equation" (Smith and Sprinks 1975)). The mild slope 


equation can be expressed in the form 


Cc 
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x and y = two orthogonal horizontal coordinate directions 
c(x,y) = the wave celerity (= o/k) 
o = angular wave frequency (defined to be 2n/T) 
k(x,y) = wave number given by the dispersion relation, 
on = gk tanh (kh) 
T = wave period 
ee) = group velocity (= 40/ak) 
o(x,y) = complex velocity potential 
g = acceleration due to gravity 
h(x,y) = still-water depth 

8. Numerical solution of this equation for the velocity potential field 
is an effective means for solving the complete wave propagation problem. The 
equation can be solved using either finite element (Berkhoff 1972 and Houston 
1981) or finite difference methods (William, Darbyshire, and Holmes 1980). 
Since transmission and reflection boundary conditions are easily implemented 
into these solution schemes, this approach is a popular one for modeling tsu- 
nami propagation and for solving problems involving the response of harbors to 
short and long waves. This method becomes computationally infeasible for 
large scale, open coast, short-wave problems because of its great expense. 
Numerical solutions of Equation 1 are only practical, as a rough rule of 
thumb, when the dimensions of the spatial area of interest are no more than 
10 times the length scale of the wave lengths being considered (Berkhoff, 
Booy, and Radder 1982). 

9. An alternative method based upon a simplification of this equation 
has recently been developed. This method alleviates the computational bur- 
den imposed by a direct solution of Equation 1. The velocity potential can 
be separated into a forward scattered and a reflected component. By neglect- 
ing the reflected part and assuming that diffractive effects in the direc- 
tion of propagation are much less than those perpendicular to the direction 
of wave advance, the following equation for the forward scattered wave can be 


derived: 
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where i = V-1 and x is now defined as the principal direction of propa- 
gation. Here, the velocity potential describes only the forward scattered 
wave field. The assumption made above, concerning the relative magnitude of 
the diffractive effects, changes the character of the governing equation from 
elliptic to parabolic. Very efficient computational techniques exist for 
solving this type of equation. Candel (1979), Radder (1979), Lozano and Liu 
(1980), Tsay and Liu (1982), Berkhoff, Booy, and Radder (1982), Booij (1981), 
and Kirby (1983) all applied this approach to study the problem of wave propa- 
gation over complex bathymetries using finite difference solution techniques. 

10. The "parabolic approximation" method, as it is called, has the fol- 
lowing disadvantage. It requires that one grid coordinate be approximately 
parallel to the predominant wave direction. This requirement can conceivably 
result in erroneous solutions to problems involving complex bathymetries where 
a dominant wave direction may not be clearly defined. Booij (1981) examined 
errors associated with the application of different parabolic approximations 
to solve the problem of oblique wave incidence over a horizontal bottom. To 
date, nothing has been documented concerning errors which may result from 
using this method to model wave incidence over arbitrary bathymetry. This di- 
rectional restriction also implies that more than one grid system may be re- 
quired in order to simulate a wide range of incident wave directions. The 
parabolic approximation method is a powerful tool for predicting linear wave 
transformations, but, it does have some deficiencies. These unaddressed prob- 
lems currently preclude its incorporation into the regional modeling system, 
as it is envisioned. 

11. The model presented in this report, RCPWAVE, is an alternative ap- 
proach for solving the open coast wave propagation problem. It addresses both 
processes, refraction and diffraction, and can be applied on a regional basis 
quite economically. The model also contains an algorithm which estimates wave 
conditions inside the surf zone. This wave breaking model is an extension of 
the work of Dally, Dean, and Dalrymple (1984) to two horizontal dimensions. 
Kirby (1983) implemented their one-dimensional breaking model into his para- 
bolic approximation model. Any short-wave propagation model integrated into 
the regional system must address the problem of wave transformation within the 


surf zone where many of the physical processes interact and move sediment. 


Wave Transformation Outside the Surf Zone 


Theoretical Basis 
12. The velocity potential function for linear, monochromatic, plane 


waves can be represented by the expression 
o = ae (3) 


where 

a(x,y) = wave amplitude function equal to gH(x,y)/2o0 

H(x,y) = wave height 

s(x,y) = wave phase function 
Here again, the velocity potential function only describes the forward scat- 
tered wave field. No considerations are given to wave reflections. By sub- 
stituting this expression for the velocity potential into Equation 1 and solv- 
ing the real and imaginary parts separately, two equations can be derived 
(Berkhoff 1976), namely, 


2 2 
: a + a + al Va-V(cec_)]}+ 1 - \vs|° =0 (4) 
ppeall Tey g 8 


v-(a-ee, Vs) = 0 (5) 


where the symbol Vv denotes the denotes the horizontal gradient operation. 

13. Together, these equations describe the combined refraction and dif- 
fraction process. Diffraction is often erroneously described as the propaga- 
tion of energy along wave crests which are defined to be perpendicular to the 
wave phase function gradient Vs . Equation 5 shows energy is still propa- 
gated in a direction perpendicular to the wave crest. Diffractiwe effects do 
change the phase function as a result of significant wave height gradients and 
curvatures. These changes cause the local wave direction to vary. If dif- 
fractive effects are neglected, Equations 4 and 5 reduce to those describing 
pure refraction in which the wave number represents the magnitude of the phase 
function gradient. 

14. Linear wave theory assumes irrotationality of the wave phase func- 


tion gradient. This property can be expressed mathematically as 


vx(Vs) = 0 (6) 


The phase function gradient can be written in vector notation as 


> 


5 
Vsil=" |Vs|cosmeham-em||Vs|msinwely (7) 
> > 
Where i and j are unit vectors in the x- and y-directions, respectively, 
and (x,y) is the local wave direction. Equations 6 and 7 can be combined 


to yield the following expression: 
Ee ( |vs| sin G), - = ( |vs| cos 0) = © (8) 


If the magnitude of the wave phase gradient is known, local wave angles can be 
calculated from Equation 8. Similarly, Equation 7 can be substituted into 


Equation 5 to yield 


es (s%c, |vs| cos 9 + a (ace, |vs| sin ) = 0 (9) 
This form of the energy equation can be solved for the wave amplitude function 
a once the wave phase characteristics Vs and 6 are known. The wave 
height can be determined and is proportional to the amplitude function, since 
wave frequency is constant. 

15. Equations 4, 8, and 9, along with the dispersion relation, describe 
the combined refraction and diffraction process for linear plane waves subject 
to the restrictions that the bottom slopes are small, wave reflections are 
negligible, and any energy losses are very small and can be neglected. The 
numerical solution scheme used to solve these equations is presented in the 
next section. These equations are assumed to be valid outside the surf zone. 
The method used to determine wave characteristics inside the surf zone is de- 
scribed later. 

Numerical solution 

16. The three governing equations (Equations 4, 8, and 9) are solved 
using numerical methods. Partial derivatives within the equations are approx- 
imated using finite difference operators. Finite difference solution methods 
require the construction of a computational grid system or mesh. Solution ac- 


curacy is directly related to resolution within the grid system. Discussion 


10 


throughout the text will refer only to grid systems comprised of constant 
sized, rectangular cells. RCPWAVE is capable of computing solutions on vari- 
ably sized, rectilinear grid systems. Technology for creating variably sized 


grids, which are compatible with the wave model, exists at CERC. 


17. Figure 1 shows nine rectangular cells which make up a small part of 


a larger mesh. Each cell has a length equal to Ax in the x-direction and 
Ay in the y-direction. The maximum values of i and j are M andy Ni % 


respectively. 


y - AXIS 


j=1TON 


i=1TOM 


x - AXIS 


Figure 1. Definition of coordinate system and grid cell 
conventions used in the model 


18. All variables which vary as a function of space are defined at the 


cell centers. For any dependent variable F , the following finite difference 


operators are used to approximate certain partial derivatives of F at the 


position (i,j): 
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Equations 11 and 13 are central differences, and Equations 10 and 12 are back- 
ward differences. All four expressions have the same order of accuracy. 
19. The magnitude of the wave phase function gradient at any point 


(i,j) is computed from the following expression, 
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This equation was derived by approximating the partial derivatives in Equa- 
tion 4 using the finite difference operators given in Equations 10 through 
13. The reason for selecting backward finite differences to approximate the 
x-derivatives, specifically the curvature, will be discussed later in this 
section. 

20. Equations 8 and 9 can both be expressed in the following general 


form: 


2 


ey LEG) (15) 


If partial derivatives in both the x- and y-directions are estimated using 
central differences about the point i-1/2,j , then an approximate form of 


Equation 15 can be written as 


Bi zi v9 [Aska itch) Het 0 i, j+1 Ltt) 20H AIG) 
The partial derivative aG/dsy at position i-1/2,j has been represented as a 
weighted average of its values at locations i,j and i-1,j . The value of 
the weighting factor W used in RCPWAVE is 1.0. This choice implies that an 
implicit solution of the equation is required. One additional approximation 
is made by using the following weighted sum: 


Pat ane +(1- 2a)F. + aF. (17) 


J iy j=! 

to represent the variable F at position i,j . Here a is another weight- 
ing parameter. The value of a is set to 0.167 in RCPWAVE. This "dissipa- 
tive interface" (Abbott 1975) is used to enhance the stability of the numeri- 
cal scheme. Substitution of Equation 17 into Equation 16 results in the 


following expression: 


F =F! 
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+ AX ( Phy +(1- wl (18) 
The finite difference formulation of Equations 8 and 9, which will be de- 
scribed next, is identical to that used by Perlin and Dean (1983). Their 
choices for W and a were 1.0 and 0.25, respectively. 
21. The finite difference form of Equation 8 is derived by substituting 


the following expressions for F and G _ into Equation 18: 


F = |vVs| sine (19) 


G = -|Vs| cos @ (20) 
These substitutions result in the following equation in which the sine of the 


local wave angle at the location i-1,j is the quantity to be determined, 


thus 


; Medien i" 
sin 6,11 5 = RS ae [elves jar sin @; 341 * (1 2a)|vs| 5 sin @; | 


+ alvs|s 54 sin eal) - HEX (les | a set COS 85 4 444 

(21) 

SS) a cos oe) - ae (lela ge cos Oriani 

-1¥S] 5 jot cos °1,3-1)] 

22. Similarly, using the substitutions 
F = acA (22) 
G = a°B (23) 
where 

A = ec, |¥s | cos 0 (24) 
B= ec, |¥s| sin 0 (25) 


the finite difference form of Equation 9 becomes 


14 


2 2 2 
a) eens (ef, serfs, er til puced ages Aaiaat oa} j-18i,5-1) 


WAX (22 (26) 


+ B - ae B 
2Ay i-1,j+1 i-1, j+1 i-1,j-1 i-1,j-1 


(1-W)Ax/ 2 2 
~ ay Suse a acapsay hy! 


This equation can be solved for the wave amplitude function, and subsequently 
the wave height, at the location i-1,j . The remainder of this section de- 
scribes the procedure used to solve the set of approximate Equations 14, 21, 
and 26. 

23. Model input (described in detail in Part IV) includes values of the 


deepwater height H direction 95 , and period T of waves to be simu- 


’ 
lated. It also ets specification of the bottom bathymetry throughout the 
grid. The wave number, which is related to the wave period and the local 
water depth through the dispersion relation, is computed at every cell. Wave 
number is used as an initial guess for the magnitude of the wave phase func- 
tion gradient. The wave celerity c and the group velocity Cy are func- 
tions of the wave period and wave number. Therefore these variables can be 
calculated at each cell. 


24. From Snell's law, 


Sime) | nee 8 (27) 
Chae meee 
fo) 
where c, is the deepwater wave celerity (defined to be gT/2n ), an estimate 


of the local wave angle is calculated everywhere. This estimate assumes that 
the bottom contours are parallel with the y-axis. If the bottom bathymetric 
contours make a known nonzero angle with the y-axis, a better first guess for 


the wave angles can be computed. The new approximation is 


sin (9 -@_) 


Q;o 
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where oF defines the contour angle. The local wave angle, deepwater wave 
angle, and contour angle follow the angle convention shown in Figure 2. The 


contour angle is an input parameter into RCPWAVE. 
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t) 


x - AXIS 0. =DEEPWATER WAVE ANGLE 


fe) 
@ =LOCAL WAVE ANGLE 
6, =OFFSHORE CONTOUR ANGLE 


Figure 2. Definition of angle conventions used in the model 


25. Wave heights at each cell are estimated as the product of the deep- 


water wave height, a shoaling coefficient rsd and a refraction coefficient 


Kee Chus 
r 
Hig= HokrXs (29) 
where 
cos 9, Vie 
“r = \eos 6 (30) 
and 


1 1/2 
(31) 


+ sinhii(ekh) tanh (kh) 


sour (: 2kh 
The dispersion relation, Snell's law, and this simple estimator of the wave 
height allow an initial guess to be made for the variables of interest 
throughout the grid system. 

26. The solution scheme implements the following marching procedure 
once initial guesses for the variables of interest have been made. Starting 
at the offshore row designated by i=M-3 , Equations 21 and 26 are used to 
compute wave angles and then heights along the entire row (from j=2 to 
j=N-1). Wave height is used interchangeably with amplitude function since 
one is directly proportional to the other. 

27. Wave angle and height solutions along a given row are solved 
iteratively because of the implicit differencing formulation used. Calcula- 
tions of the wave angle (actually the sine of the wave angle) and the wave am- 
plitude function are reiterated until the average change (along a row) in each 
variable from one iteration to the next is less than some tolerance. These 
convergence criteria, 0.0005 for wave sines and 0.001 ft* (or a metric equiva- 
lent) for wave heights, are suggested values for prototype applications. They 
can be easily changed by modifying the source code using the method outlined 
in Part IV. 

28. This solution considers only refraction since the wave number k 
is used as an estimate of the magnitude of the phase function gradient. Equa- 
tion 14 is then used to compute the true magnitude of the wave phase gradient. 
This new wave number accounts for the effects of diffraction. Backwards dif- 
ferences are used to approximate the x-derivatives because they require only 
information which has already been computed. Next, Equations 21 and 26 are 
again solved in order to compute the wave angles and heights using these new 
wave numbers. This procedure is repeated along the row under consideration 
until the change in new wave number, from one iteration to the next, is less 
than 0.5 percent of the newly computed value. This condition must be met at 
each cell along the row. As a row of new wave numbers is computed, the values 
are filtered in the y-direction using the method of Sheng, Segur, and Lewellen 


EEE 


* To convert feet to metres use a conversion factor of 0.3048. 
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(1978). This filter removes cell-to-cell oscillations introduced as a result 
of the differencing scheme used to compute the new wave numbers. Row-by-row 
marching proceeds until solutions are computed along row i=2. 

29. Lateral boundary conditions for a row are specified at the conclu- 
sion of calculations for that row. The value of all variables at cells j=N 
and j=1 are set equal to their values at cells j=N-1 and j=2 , respec- 
tively. This boundary condition implies that the change in the variable in 
the y-direction is zero. The condition is most valid when the bathymetric 
contours are nearly straight and parallel to the y-axis. For this reason it 
is recommended that users orient their grid system so that the y-axis is 
nearly parallel to bottom contours along the lateral boundaries. 

30. Boundary conditions along the seaward extent of the grid are used 
to initiate the shoreward marching algorithm. They are computed from deep- 
water wave input supplied by the user along with the following assumption. 
Bottom contours extending from the offshore grid row i=M out to deep water 
are assumed to be straight and parallel to a line making an angle of 8, with 
the y-axis. In other words, Snell's law is assumed to be valid from deep 
water to the outer boundary of the grid system. No inshore boundary condi- 
tions (along row i=1) are required because of the foreward marching solution 


scheme. 
Wave Transformation Inside the Surf Zone 


Theoretical basis 

31. Waves approaching the very nearshore zone tend to steepen and even- 
tually break because of decreasing water depths. Shoreward of this breaking 
point dissipative energy losses due to turbulence strongly influence the wave 
height. Linear theory does not allow for prediction of the breaker location 
nor for wave transformation across the surf zone. Instead, empirical and ap- 
proximate methods must be used to describe the breaking process. 

32. The first aspect to consider in surf zone transformation of waves 
is incipient wave breaking. Iwata and Sawaragi (1982) reviewed many criteria 
for determining wave characteristics at the breaking point. One which is ap- 
pealing because of its basis on wave physics defines breaking as the point 
when the particle velocity at the wave crest exceeds the wave celerity. The 


following formulas: 
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and 


b L 


enh, 
H. = 0. 142-0 ‘tanh ( (33) 
b b 


H, = breaking wave height 
= water depth at breaking 

Ly = wave length at breaking 
by McCowan (1891) and Miche (1944), respectively, are based on this criterion. 
Equations 32 and 33 were derived for solitary and periodic shallow-water 
waves, respectively, in water of uniform depth (Iwata and Sawaragi 1982). A 
breaking height predictor based on wave energy flux was developed by Komar and 


Gaughan (1972) and is given by 


H. = (gy (mre) (34) 


where «* is a dimensional coefficient equal to 0.39. Field and laboratory 
data have shown this predictor to be quite accurate. Other incipient breaking 
criteria have been developed by fitting empirical relationships to breaking 
wave data. These methods are not derived from any theoretical considerations 
of wave physics, yet results derived using them agree very well with observed 
data. The most widely used of these criteria are those developed in Equa- 
tions 35, 36, and 37 by Le Mehaute and Koh (1967), Goda (1970), and Weggel 


(1972), respectively. These criteria are as follows: 


H -1/4 
% s10) WA 
H, = 0.76 HL| 7 m (35) 
fo) 
where 
Lo = deepwater wave length 
m = bottom slope 
hy, 4/3 
H, = 0.17L5)1 - exp) -1.5 --(1 + 15(m)) (36) 
fe) Ly 
and 
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H, = ——2.— (37) 


where 


o(-19m)] 
(-19.5m)] 


43.7511 
1.56/11 +e 


a 


b 

33. All of these empirical methods give reasonable approximations of 
the incipient breaking wave height. In choosing a criterion for inclusion 
into RCPWAVE, the following factors were considered. The criterion should ac- 
count for bottom slope and wave period since field and laboratory tests show 
these parameters to be important (Iwata and Sawaragi 1982). It should not 
depend on deepwater parameters alone because the transformation algorithm must 
be capable of modeling multiple breaking and reformation. The Weggel (1972) 
and Goda (1970) criteria satisfy both requirements. The former was selected 
for inclusion into RCPWAVE. 

34. Once the incipient breaking point is defined, a mechanism is needed 
to transform the breaking wave across the surf zone. Historically the wave 
height has been assumed to be proportional to the local water depth throughout 
the surf zone. The constant of proportionality was assumed to be about 0.8. 
Field and laboratory data have shown that this approximation consistently 
overestimates actual wave heights within the surf zone (Dally 1980 and Thorn- 
ton and Guza 1982). Other formulations have been developed and applied suc- 
cessfully by researchers to calculate surf zone wave heights. Most are of the 


following form: 


9(Ec_) 
ae a sas (38) 


ax = 


where Ec, is the energy flux associated with the breaking wave and 6 isa 
term representing the rate of energy lost due to bottom friction, turbulence, 
and other dissipative processes. 

35. Divoky, Le Mehaute, and Lin (1970) considered the dissipation in a 
breaking wave to resemble that in a hydraulic jump. They assumed that the 


change in energy flux could be approximated by 
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a(Ec_) (y, =Y4 \3 
= 89 | hy ¥ (39) 


where 
p = water density 

Q = flow across the jump 

Y, = water depth on the high end of the jump 

Y5 = water depth on the low end of the jump 
For water of uniform depth, they assumed that Y5 =whi+! Ale vand Yy = h + BH 
where 8 is a coefficient related to the percentage of the wave height 
covered by foam. Battjes and Janssen (1978) also used a hydraulic jump repre- 
sentation of wave breaking. They used the following expression: 


9(Ec_) ye 
ax =a* 0210 To (40) 


where a* is a coefficient of order one. Mizuguchi (1980) modeled surf zone 


energy loss as, 


a(Ec_) 2 
g kH 
ax = GRE (4) ) 


where Vis is a coefficient of turbulent eddy viscosity. Results derived from 
this model compared quite well with experimental data, but any physical rea- 
soning behind the use of the eddy viscosity formulation is, as Mizuguchi 
(1980) states, obscure. Horikawa and Kuo (1966) developed an analytical 
scheme using second order solitary wave theory and theoretical expressions for 
dissipation due to bottom friction and turbulence. Their governing equation 
contained two coefficients, one for each dissipative mechanism. 

36. The transformation algorithm selected for use in RCPWAVE (Dally, 
Dean, and Dalrymple 1984) uses the same energy flux basis as the models men- 
tioned above. However, instead of using a hydraulic jump to represent energy 
loss in a single breaking wave, the form of the hydraulic jump energy loss is 
used to approximate losses across the entire surf zone. Through analogy with 


energy flux in a channel, the following equation is postulated: 


a(Ec_) 


pl oh ape, 
Es lec, ; (F,),| (42) 


fall 


where 


k = rate of energy dissipation coefficient (set equal 
to 0.2 in RCPWAVE) 


(cg) = stable level of energy flux that the transformation 
S process seeks to attain 


The right-hand side of Equation 42 is simply a dissipation term. The sub- 
script s is used to denote the stable level of some quantity. Substituting 
the linear wave theory estimate for E (E = 0.125 ae) into Equation 42 re- 


sults in the following expression: 


9(H@c_) ee 
ee =i! Cg i (43) 


37. Various field (Thornton and Guza 1982) and laboratory (Horikawa 
and Kuo 1966) experiments have shown that, well into the surf zone, the wave 
height tends toward a stable value which is proportional to the local water 


depth. This relationship can be expressed as 


where 


= 
! 


= stable wave height 
Y = proportionality coefficient (set equal to 0.4 in RCPWAVE) 


Equation 43 can now be rewritten as 
2 
a(H © ) 
ax h g - 


38. This surf zone wave transformation model can be incorporated into 
the conservation of wave energy equation (Equation 5) by simply adding the 
dissipation term D to the right-hand side. The function D must now repre- 
sent dissipation in the direction of wave propagation. Also for dimensional 
consistency, the term D must be multiplied by the wave celerity and the mag- 
nitude of the wave phase gradient, and the wave height must be replaced by the 
wave amplitude function. In vector notation, the energy equation becomes 


v-(2°c0,95) = = | sPespiv9 - (Ga) haces \vs| é (46) 
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This equation can be thought of as being valid both inside and outside the 
surf zone. Outside, the coefficient « is zero, and the equation reduces to 
Equation 5. 

39. Discussion relating to wave transformation within the surf zone has 
addressed the problem of determining wave heights. The problem of wave phase 
must be addressed also. Diffraction effects are assumed to be negligible 
inside the surf zone. Therefore, the wave number k is assumed to accurately 
represent the magnitude of the wave phase function gradient. The linear wave 
theory assumption that the waves are irrotational also will be assumed to re- 
main valid inside the surf zone. Consequently, wave angles inside the surf 
zone are computed in the same manner that was used outside the surf zone. 
Numerical solution 

40. The numerical procedure for computing wave angles inside and out- 
side the surf zone is the same. This section documents only the solution 
scheme used to determine breaking wave heights. The finite difference form of 
the wave energy equation outside the surf zone (Equation 26) can be expressed 


in the following form: 


ac ty F + Ax G (47) 
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With the inclusion of the dissipative term, Equation 47 becomes 
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a. a = ea a ee (48) 
where D* represents the finite difference form of the dissipation term on 
the right-hand side of Equation 46. Reiterating, the dissipation term is an 
average value along the wave path. The wave path is determined by the local 
wave angle at the position i-1,j which has already been computed. There- 
fore, the average along the path is an average of information at cell i-1,j 
and another cell whose position is denoted by ikey,jkey . The procedure used 
for determining the location of this cell will be presented later. 

41. The term D can be written in finite difference form as 


|vs| + (ace | Vs | 
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where 


With some algebra, Equation 48 can be reorganized so that the amplitude func- 
tion at the position i-1,j only appears on the left-hand side of the equa- 
tion. Therefore, the energy equation inside the surf zone can be numerically 
solved using the same procedure which was used to solve it outside the surf 
zone. 

42, The location of the cell denoted ikey,jkey is found using the 
following procedure. "Areas of influence" are determined by extending lines 
from the center of the cell i-1,j to the midpoints between the surrounding 
cell centers (Figure 3). Angles are computed from the x-axis to these radial 
lines. The local wave angle calculated at cell i-1,j is compared to each of 
these angles in order to determine the nearest, prior cell along the wave 
path. For example (refer to Figure 3), if the local wave angle is greater 
than 6 but less than 6 


2 
ikey = i and jkey = j+1. 


1? then cell i,j+1 is the cell of influence and 


43, A flow chart describing the wave height computation is shown in 
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CELL 
CENTERS 


Figure 3. Schematic drawing showing cell 
of influence conventions used in the wave 
breaking scheme 


Figure 4. The wave amplitude function is computed from the energy equation 
assuming no dissipation. The amplitude function is converted to wave height 
and compared to the stable wave height yh . If the wave is less than or 
equal to this stable level, the wave is either inside the surf zone, having 
been transformed to a state below the stable level, or it is outside the surf 
zone. In either case, no further transformation is needed. If the wave height 
is greater than yh , then additional wave transformation may occur. The cell 
of influence is located and tested to determine whether or not the wave has 
experienced prior breaking. If the wave is undergoing transformation in the 
cell of influence, it continues to be transformed. If the wave in the cell of 
influence is not being transformed, the local wave height is checked against 
the incipient breaking height criterion. If the height exceeds the allowable 
value, wave dissipation begins. The accuracy of the surf zone wave transfor- 
mation model has been verified using laboratory data of Horikawa and Kuo 
(1966) and Izumiya (1984). Comparisons between model results and these data 


are described in Part III. 
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Figure 4. Flow chart of the wave breaking scheme 
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PART III: MODEL VERIFICATION 


General Comments 


44, Comparisons between model results and observed data were used to 
verify the model. Both laboratory and field data were used in these tests. 
The ability of RCPWAVE to simulate wave transformation outside the surf zone 
was checked using data collected during a laboratory experiment conducted by 
Berkhoff, Booy, and Radder (1982) and using prototype data obtained during a 
field experiment at the CERC Field Research Facility (FRF) in Duck, North 
Carolina. The next two sections describe these comparisons in detail. Fig- 
ures pertaining to the laboratory and field verification are located in Appen- 
dixes A and B, respectively. Only laboratory data were used to verify the 
surf zone wave transformation part of the model. These data were collected 
during one-dimensional flume tests performed by Horikawa and Kuo (1966) and 
Izumiya (1984). Both experiments considered only breaking of monochromatic, 
plane waves. The former experiment investigated wave transformation on a 
plane beach only; the latter involved tests using plane, stepped, and barred 
beaches. These comparisons are discussed in the last section. Figures per- 


taining to the wave breaking verification are located in Appendix C. 


Elliptical Shoal Case: Comparison with Laboratory Data 


45, Berkhoff, Booy, and Radder (1982) performed a wave tank experiment 
in which wave conditions resulting from the propagation of a monochromatic, 
plane wave train over complex bathymetry were measured. The bathymetry con- 
sisted of an elliptical shoal superimposed on a plane beach. This geometry is 
shown in Figure A1. The shoal acts as a lens and focuses incoming wave energy 
into a strong convergence zone. This experiment provided a set of data which 
could be used to verify the ability of RCPWAVE to compute accurate solutions 
to problems where refraction theory breaks down. Wave height data were col- 
lected at many locations within this zone along the eight cross sections shown 
in Figure Al. RCPWAVE was used to simulate the same experiment. Simulated 
wave heights were then compared with observed data. 

46. A finite difference grid mesh, which measured 25 m in the 


x-direction and 20 m in the y-direction was constructed. Rectangular grid 


eth 


cells measuring 0.25 and 0.20 m in the x- and y-directions, respectively, were 
used for the simulation. The center of the shoal was located at the point 
with x- and y-coordinates of 15 and 10 m, respectively. The equation de- 


scribing the outer extent of the shoal is given by 


y" 2 oA 2 

Saline) ao (50) 
where x' and y' denote a local coordinate system whose origin is located 
at the center of the shoal. Also this coordinate system is rotated 20 deg in 
the clockwise direction relative to the x- and y-coordinate system. All length 


scales for both coordinate systems are measured in metres. Water depths, de- 


noted by h and measured in metres, were calculated from the following three 


equations, 
hig=s 0-4 Sya ora yy ii5)04 (51) 
h = 0.45 - 0.02 (5.84 - y')(outside the shoal boundary) (52) 


2 2 1/2 
he =0.45) = 0.0205. 845 — y's) 45053) =—9 025i = ‘2 - (e (53) 
(inside the shoal boundary) 


An incident wave with a period of 1.00 sec, a height of 1.06 cm, and a direc- 
tion of approach parallel to the x-axis was used as the deepwater boundary 
condition for the simulation. 

47. Figures A2 through A4 show comparisons between the experimental 
data (open circles) and model results (solid line) for all eight profiles. 
Results are presented as profiles of relative wave height (observed wave 
height divided by incident wave height). A refraction analysis (presented in 
Berkhoff, Booy, and Radder (1982)) shows that a caustic occurs at a point be- 
tween profiles 2 and 3. Results show that the model is quite capable of accu- 
rately simulating diffractive effects at, and beyond, the caustic location. 
Consequently, it is a much more powerful tool for simulating wave transforma- 
tion than refraction theory alone. Berkhoff, Booy, and Radder (1982) and 
Kirby (1983) also numerically simulated this experiment in order to test their 
parabolic approximation models. Accuracy of their model results, using linear 


wave theory, is comparable to that obtained using RCPWAVE. Results obtained 


by Kirby (1983) are shown also in the figures (as a dashed line). Kirby 
showed that much of the discrepancy between simulated and observed data could 
be eliminated if a nonlinear wave theory were incorporated into the model. He 
incorporated Stokes' second order theory (which he showed was valid for this 
experiment) into his model, and the results also showed that nonlinear effects 
became increasingly important after the waves pass profile 3. 

48. An interesting aspect of the RCPWAVE results is evident in pro- 
files 3, 4, and 5. The lobed features in the wave height variation are 
smoothed by the model. The cause of this smoothing is not known. It may be 
caused by the dissipative interface or the point-to-point filter used in the 
numerical scheme. The side lobes seem to be related to the occurrence of an 
amphidromic point where the wave phase becomes multivalued and the wave height 
variation contains a discontinuity. The solution scheme forces the magnitude 
of the phase function gradient to be single valued. This may also cause the 
local smoothing. The model is intended for use in open coast, prototype ap- 
plications. For these types of problems, this smoothing property of the model 


can certainly be tolerated. 


CERC's FRF Cases: Comparisons with Prototype Data 


49. In addition to the laboratory verification, RCPWAVE also was veri- 
fied using field data collected at CERC's FRF in Duck, North Carolina. Bot- 
tom bathymetric contours in the area are generally straight and parallel to 
the coastline except in the immediate vicinity of the research pier. The 
pier's presence has caused the formation of a deep scour hole along much of 
its length. The complicated bathymetry, which has resulted from this hole, 
was one reason for selecting the FRF for field verification. Hubertz (1981) 
showed that a ray model using refraction theory alone proved incapable of sim- 
ulating observed conditions. Hubertz (1982) also showed that a short wave 
model which includes diffractive effects in its governing equations (the Sys- 
tem 21 Mark 8 proprietary model developed by the Danish Hydraulic Institute) 
could accurately simulate wave propagation in the vicinity of the pier up to 
the breaking point. 

50. Another reason for selecting the FRF was the availability of 
wave data, both offshore and along the pier. During October 1982, an exten- 


Sive, 1-month field data collection program was undertaken. Two storms 
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occurred during the month, one lasting from 10-13 October and the other from 
23-26 October. The latter was the more severe event. Prototype data col- 
lected during this month included bathymetric surveys conducted on October 16 
and 27, continuous and synchronous wave data from six Baylor gages along the 
pier and an offshore Waverider buoy, radar imagery of the sea surface, and 
continuous tide elevation data. All wave data were available in a spectral 
form (energy density as a function of frequency). Estimates of significant 
height and peak period were provided also. Wave directions could be deter- 
mined from the radar imagery. The availability of clear radar pictures was 
ultimately the discriminating factor in selecting which time periods were con- 
sidered for verification purposes. 

51. Six cases, or time periods, were ultimately chosen using the fol- 
lowing process. As stated above, only those times were considered when clear 
radar imagery was available. The number of candidate times were further re- 
duced by requiring that wave data be available within approximately 1 hr of 
the time the image was taken. Next, the wave spectum from the offshore gage 
was checked for spectral shape. Only times with single-peaked, fairly narrow 
banded spectra were considered. Six cases which met all the criteria are 


shown in Table 1. 


Table 1 
Summary of Field Data Used to Verify RCPWAVE 


Time of 
Test 9 Radar Bathymetry 
Case Time Hy AB fe) Tide Imagery Survey 
Number Date G.m.t.* m sec deg m G.m.t. Date 
1 10-13-82 1300 1.95 13.21 -23 0.12 1115 10-16-82 
2 10-13-82 1400 We8i7, 14.22 -25 -0.18 VAS 10-16-82 
3 10-15-82 1210 0.78 12.34 -25 0.78 1130 10-16-82 
4 10-17-82 1200 1.56 6.87 43 0.91 1120 10-16-82 
5 10-25-82 1900 3.10 12.34 -25 0.68 1950 10-27-82 
6 10-25-82 2000 2.95 12.34 -25 0.63 1950 10-27-82 


* G.m.t. denotes Greenwich mean time. 


52. For each case the date and time of the recorded wave data are 
shown. The parameters Hy vel and oF are the deepwater wave characteris- 
tics used as input into RCPWAVE. The value of T was chosen to be the best 
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estimate of the peak spectral period and was made using data from the Wave- 
rider buoy and the seaward-most Baylor gage located at the end of the pier. 
The deepwater wave angle was computed using this peak period and a wave direc- 
tion estimated from the radar imagery. The time at which the radar image was 
taken also is given. Inherent in the procedure used to calculate the wave 
parameters in deep water is the assumption that the bottom contours seaward of 
the pier are straight and parallel. This assumption is quite reasonable for 
this stretch of coastline. The deepwater wave height represents a significant 
height and was computed using the peak period, deepwater angle, and the sig- 
nificant height recorded by the offshore Waverider buoy. The recorded wave 
spectra at the offshore gage are shown for all six cases in Figure B1. 

53. Bottom bathymetry is required as model input. The total water 
depth matrix, used in the model, is computed by simply adding some tidal ele- 
vation to each depth value. Depth values were taken from one of two surveys 
shown in Figure B2. The particular survey used for each verification case is 
shown in Table 1. The tidal elevation (relative to mean sea level (MSL)) also 
is given in Table 1. The areal extent of each survey is identical, covering 
1,200 m in the y-direction and 900 m in the x-direction. The orientation of 
the survey axes was adopted for use in constructing the model grid system. 

The x-axis is parallel to the FRF pier. Actual depth values (relative to MSL) 
were provided for each cell of a grid comprised of 75 cells in the x-direction 
and 50 cells in the y-direction. This grid completely encompasses the sur- 
veyed region. Cell dimensions are 12 and 24 m in the x- and y-directions, 
respectively. 

54. Comparisons between simulated and observed significant wave heights 
along the pier are shown for each case in Figures B3 and BY. For these tests, 
the model is being used to propagate some amount of energy (here, designated 
by the significant height) with a single frequency and some mean direction. 

By requiring that the radar imagery be clear and contain only a unidirectional 
wave train, only waves which are nearly planar (long-crested) are being con- 

sidered. Since most of the wave spectra are narrow-banded (with the exception 
of Case 4), the cases being considered represent nearly monochromatic condi- 

tions. Therefore, assumptions inherent in the model's governing equations are 
essentially upheld, and the model should be able to simulate these conditions. 
Results show that RCPWAVE accurately predicts wave propagation for these types 


of wave conditions over a complex bottom. In all cases the scour hole causes 
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wave energy to be propagated away from the hole causing a reduction in wave 
height along the pier. Results from Hubertz (1982) substantiate this observa- 
tion. The trend of the wave height variation along the pier for Case 4 also 
is accurately simulated. The magnitude of the simulated significant height 
consistently underestimates observed data by about 0.1 m. The closeness of 
fit for this case, which involves a wider spectral bandwidth, suggests that 
the model may provide a useful method for predicting spectral wave transforma- 
tions if the waves have a small directional spread. More research is needed 


to test this hypothesis. 


Verification of the Wave Breaking Scheme: Comparisons 
with Laboratory Data 


55. A number of numerical simulations were performed in order to verify 
the capability of RCPWAVE to predict wave breaking and surf zone wave trans- 
formation. Model results were compared with data from laboratory experiments. 
Two sources of data were used. The first source was the original work of 
Horikawa and Kuo (1966) along with additional information concerning that work 
found in Dally (1980) and Dally, Dean, and Dalrymple (1984). Horikawa and Kuo 
studied the transformation of waves in the surf zone using two wave tanks and 
four different uniform bottom slopes. Two slopes, 1:20 and 1:30, were tested 
in a wave tank 17 m long and 0.6 m deep. The remaining slopes, 1:65 and 1:80, 
were tested using a larger tank which measured 75 m in length and 1.2 m in 
depth. The second source of experimental data was the dissertation work of 
Izumiya (1984), who used a smaller wave tank for his tests. All water depths 
considered in the experiments were less than 0.3 m, and all breaking wave data 
were collected within 7 m of the dry "beach." He investigated wave transfor- 
mation over three different bottom configurations: a plane beach, a stepped 
beach, and a barred beach. All slopes were 1:20, including the back slope on 
the barred beach. 

56. Wave and bathymetry conditions for each laboratory test used in the 
model verification process are shown in Table 2. The table shows the data 
source, an arbitrarily assigned test case number, the type of beach consid- 
ered, and the beach slope. The following wave parameters for each test are 
given also: (a) the wave period, (b) the deepwater wave height (if it were 
available), (c) the incipient breaking wave height, and (d) the water depth at 


breaking. 
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Table 2 
Summary of Laboratory Data Used to Verify the Wave Breaking Scheme 


at Hy Ay hy 
Test Case Bathymetry Slope sec cm _em_ mem 

Horikawa and Kuo 1 Plane beach 1:20 1.40 -- 13.2 5c 
Horikawa and Kuo 2 Plane beach 1:30 2.20 -- 9.1 10. 
Horikawa and Kuo 3 Plane beach 1:30 2.20 -- 12.2 16 
Horikawa and Kuo 4 Plane beach 165 56 14.4 Wee 19. 
Horikawa and Kuo 5 Plane beach 165 1.56 24.7 2 teil 35). 
Horikawa and Kuo 6 Plane beach 1:80 1.40 -- leet 25 
Izumiya 7 Plane beach 1:20 1.19 6.0 Val 9 
Izumiya 8 Stepped beach 1:20 0.96 7.9 8.5 9 
Izumiya 9 Barred beach 1:20 0.95 6a3 6.9 


57. Two series of verification tests were conducted. The first series 
involved model simulations which were initiated at the break point determined 
from the data in Table 2. The purpose of these tests was to check the valid- 
ity of the surf zone transformation part of the breaking scheme. Figures Cl 
through C6 show plots comparing simulated results with Horikawa and Kuo's 
data. The horizontal and vertical scales vary from plot to plot. Model re- 
sults generally show good agreement with observed data. The model tends to 
underpredict the wave heights very close to shore because the model does not 
consider the effects of wave setup. Setup would increase the total water 
depth here and allow for larger simulated wave heights (Dally 1980). Horikaw 
and Kuo Test 2, having data points 1.0 and 0.5 m from shore, clearly illus- 
trate this model deficiency. However, neglecting wave setup has little effec 
on model accuracy in the remainder of the surf zone. Steep slopes, 1:20 and 
1:30, produce a much more rapid decrease in wave height across the breaking 
zone than do the milder slopes, 1:65 and 1:80. Qualitatively, this is what i 
usually observed in the field. Plunging waves are more likely to occur on 
steep slopes and spilling breakers on mild slopes. Results suggest that the 
breaking model is applicable under both types of breaking conditions. 

58. Figures C7 through C9 show comparisons between model results and 
Izumiya's data, again for surf zone transformation only. Test 7, a plane 


beach case, shows very good agreement. The wave height variation across the 
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surf zone is characterized by a strong gradient, typical of steep slopes. 

Test 8 was conducted using a stepped beach. The model reproduces the same 
general wave height variation exhibited by the data but not as accurately as 
in the plane beach cases. The data show a more rapid decay immediately inside 
the break point, and the decay seems to be toward a stable height less than 
that predicted by the model. Some oscillations appear in the data, possibly 
from reflections, which make an estimation of the stable height difficult. 

The values of « and y used in the model's breaking scheme were selected to 
provide a "best" fit to both data sets. Better agreement could have been ob- 
tained for individual data sets by using some other values. There is also 
evidence that the values of « and y are dependent upon the beach slope 
(Dally, Dean, and Dalrymple 1984). For simplicity, a constant value for both 
parameters is used in the model. The final comparison is for the barred beach 
case. Again, good agreement is obtained for the overall shape of the wave 
height decay. The model miscalculates the location of the second break point 
by about 0.3 m, but the simulated decay after this point closely parallels the 
data. 

59. A second series of comparisons, using the same laboratory data, 
were made. Results from these comparisons are provided to show the capability 
of the model to simulate the entire shoaling and breaking process, including a 
determination of the break point. Previously, only the accuracy of the decay 
mechanism was examined. RCPWAVE was used to simulate wave propagation from 
generation through to breaking. These simulations were performed for those 
cases in which the deepwater wave height was given by the authors or could be 
computed. The deepwater wave height is used as model input. Figure C10 shows 
the model results from Horikawa and Kuo Test 5. The model could not shoal the 
wave up to the observed incipient breaking height; however, the simulated surf 
zone transformation matches well with observed data. 

60. Results from similar comparisons between model results and labora- 
tory data for the Izumiya Tests 7, 8, and 9 are shown in Figures C11 through 
C13. In all cases the model failed to reproduce the wave shoaling peak imme- 
diately prior to breaking. The model does decay the wave correctly for Tests 
7 and 8. In Test 9 the model misses the first break point entirely. This 
case was included to illustrate the dependency of accurate surf zone simula- 
tions on the validity of the incipient breaking criterion used. 


61. For the monochromatic wave conditions considered in this 
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verification, and probably narrow-banded spectra as well, the model is inca- 
pable of predicting the full extent of the shoaling peak prior to breaking, 
probably due to the use of linear wave theory. This limitation, coupled with 
any errors in the incipient breaking criterion, may cause the location of the 
break point to be erroneously specified (e.g., Izumiya Test 9). However, the 
surf zone transformation part of the model seems to be very accurate for plane 
beaches and quite reasonable for stepped and barred bathymetry. Consequently, 
the model is usually capable of simulating the general form of the wave height 
variation across the breaker zone even though there may be a slight horizontal 
shift between observed and simulated distributions. This is substantiated by 
the results presented above. Inclusion of this breaking scheme into the 
model, despite its shortcomings, is an obvious improvement over methods which 
specify the surf zone wave height to be proportional to the local water depth 
everywhere, especially for complex bathymetry involving stepped or barred 
beaches. The accuracy of this model for use in predicting surf zone charac- 
teristics under wide-banded spectral wave conditions is questionable. The 
model assumes wave breaking begins at a well defined point. This assumption 
is most valid under extremely narrow-banded spectral conditions but is vio- 


lated otherwise. 
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PART IV: MODEL EXECUTION ON THE CONTROL DATA 
CORPORATION COMPUTING SYSTEM 


General Comments 


62. The RCPWAVE model can be run in a remote batch mode on the Control 
Data Corporation (CDC) computing system. Batch jobs comprised of job control 
language (JCL) allow users to instruct the computer to perform a series of 
tasks including file manipulation, compilation, and execution. RCPWAVE can be 
run on two CDC computers, the CYBER 865 and CYBER 205. Costs are less expen- 
sive on the CYBER 865 machine (cost comparisons will be given in Part V); how- 
ever, available memory within this computer limits its use to applications in- 
volving fewer than 9,025 grid cells. The CYBER 205 has more available memory, 
but it is more expensive to use. It is recommended that users restrict their 
use of the CYBER 205 to model applications where the total number of grid 
cells exceeds 9,025. 

63. JCL used to run the wave model on the CYBER 865 differs from that 
used to run it on the CYBER 205. The next two sections explain CYBER 865 JCL 
and CYBER 205 JCL, respectively. Then the following procedures are explained: 
submitting batch jobs on both machines, checking job status, and accessing job 
output. The last two sections give a general description of necessary input 
files and all output files which are generated in the course of job execution. 
Part V of the text covers examples of actual applications of RCPWAVE run on 
both CDC machines, and it includes specific examples of the aforementioned 
files. 

64. It is essential for the user to be familiar with XEDIT* or another 
text editor available on the CDC system. This knowledge is required in order 
to make any changes to the JCL files or to examine an output file in a text 
editor mode. The text editor also is required to create input data files. A 
limited knowledge of the "UPDATE"* method for making changes to the model 
source code also is recommended. Discussion concerning modifications to the 
source code and examples documenting this procedure are given in subsequent 


sections. 


* Cybernet services manuals are available from CDC titled "XEDIT, Extended 
Text Editor," Manual No. 76071000C, and "UPDATE, Version 1 Reference 
Manual," Manual No. 60449900B. 
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CYBER 865 Job Control Language 


65. An example of JCL used to run RCPWAVE on the CYBER 865 is given in 
Figure 5. The file containing this JCL can be obtained by logging into the 
CDC computing system and typing 


GET , RCP865J/UN=CER@Q2 


This action will create a local file called RCP865J on the user's work space. 


In order to save the file permanently, type 


SAVE, RCP865J 


Each line or "card image" of RCP865J contains a specific instruction to the 
computer. Note the brief description of important input and output files un- 
der the "COMMENT" section. These files will be discussed in greater detail 
later in Part IV, and specific examples will be given in Part V. At this 
point, it is only important to know what each file contains. A line-by-line 


description of the JCL in RCP865J is given in Appendix D. 


Submitting the CYBER 865 Batch Job 


66. A series of commands (JCL) has been assembled into a batch job that 
instructs the computer to gather input, compile and execute the wave model, 
and save output. The computer is actually given this batch of commands by 


using the following form of the SUBMIT command: 
SUBMIT,<JCL FILENAME> ,T 
The <JCL filename> of the file shown in Figure 5 is RCP865J; therefore the 
following command is given: 
SUBMIT, RCP865J,T 
The parameter T tells the computer to return certain output to the user's 


terminal once the job is completed. When a batch job is submitted, the com- 


puter responds by giving the job an identifier similar to the one shown below: 


15.21.05. SUBMIT COMPLETE. JOBNAME IS AKQBRMA 
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The job identifier (or job name) is denoted by a sequence of seven letters or 
numbers, here AKQBRMA. The first four characters are uniquely associated with 
a particular user number. The last three characters identify the job which 
was submitted. 

67. In order to check the progress of a batch job submitted to the 


CYBER 865, the user can issue the following command: 


ENQUIRE, JN=RMA 


The computer will respond in one of the following ways: 


a. AKQBRMA IN INPUT QUEUE--This means that the batch job AKQBRMA 
is in the input queue; i.e., the job has not begun execution. 


b. AKQBRMA EXECUTING--This means that the batch job AKQBBQF is 
executing; i.e., the computer is processing the JCL commands. 


c. AKQBRMA IN PRINT - TERMINAL QUEUE--This means that the job 
AKQBRMA is in the PRINT-TERMINAL queue; i.e., the batch job has 
been completed, and the output is available to the user. 


68. By using the T parameter on the SUBMIT command, certain infor- 
mation concerning the batch job is returned to the terminal upon its comple- 
tion. This information includes, in the following order, a copy of the diary 
of events, actual user cost information, CDC scheduling information, a com- 
piled listing of RCPWAVE, and abnormal execution information. The following 


QGET command is used to look at this information: 
QGET, XYZ 
The user must specify only the last three letters of the computer-generated 


job name (xyz). For example, if the batch job AKQBRMA is in the PRINT- 
TERMINAL QUEUE, the following command would be entered: 


QGET, RMA 


The computer would respond with 


AKQBRMA ATTACHED 


Then, the text editor could be used to look at this output by typing 


XEDIT, AKQBRMA 


69. The same information accessed with the QGET command resides on the 


permanent files RCPDAYF and RCPOTPT. These permanent files can be viewed in 
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three different ways. They can be accessed with a GET command and then exam- 


ined via the text editor with the following information: 
GET , RCPDAYF 
XEDIT, RCPDAYF 
and 
GET , RCPOTPT 
XEDIT, RCPOTPT 


Alternatively, the files can be listed at the user's terminal by typing the 


following commands: 
OLD, RCPDAYF 
LIST 
and 
OLD, RCPOTPT 
EST 
The files can be listed also at a remote printer site by entering 
ROUTE , RCPDAYF , DC=PR, UN=SC 
and 


ROUTE , RCPOTPT , DC=PR, UN=SC 


where SC is the site code of the device which is to receive the printed file. 
The printed output file RCPPRNT and the input files RCPDATA, RCPUPDT, and 
RCPDEPT also can be listed in any of these three ways. 


CYBER 205 Job Control Language 


70. Running RCPWAVE on the CYBER 205 is both more complicated and ex- 
pensive than on the CYBER 865. Before work can be done on the CYBER 205 a 
separate CYBER 205 user number and password must be obtained. Additionally, a 
direct access file must be created on the user's account which contains vali- 
dation information (USER and CHARGE cards for his/her CYBER 205, and CYBER 865 
accounts). This file will be called AF205. 


ho 


71. An example of JCL which creates this file is shown in Figure 6. 
The file containing this JCL can be copied onto the user's local file space 


and then saved permanently by typing in the following commands: 
GET , AFJCL/UN=CER@Q2 


SAVE, AFJCL 


SENDJOB. 
USER(U=(205 USER NAME) ,PA=(205 PASSWORD) )ADY 
RESOURCE ( JCAT=P3,TL=200) 

CHARGE( (CHARGE NO.), PROJECT NAMED) 

TV, 10+. 

PURGE (AF205) 

TV, 4. 

COPY, INPUT ,AF205. 

DEFINE ,AF205. 

USER(<B65 USER), <B65 PASSHORD) ,KOE) 

CHARGE( <CHARGE NO.), (PROJECT NAMED) 


Figure 6. File AFJCL, JCL to 
create validation information 
for CYBER 205 usage 


By submitting this JCL file for execution, the file AF205 will be created. 
Certain information must be entered into this file before it is submitted: 

the user's CYBER 205 user number and password (<205 USERNAME> and <205 
PASSWORD>, respectively), the user's charge information (<CHARGE NUMBER> and 
<PROJECT NAME>) with which he/she logged in, and the user's CYBER 865 user 
number and password (<865 USERNAME> and <865 PASSWORD>, respectively). Once 
these changes have been made, the JCL file can be submitted with the following 


command: 
SUBMIT, AFJCL,T 


Successful completion of this batch job creates the file AF205 which is re- 


quired in some of the CYBER 205 commands. 
72. An example of JCL used to run RCPWAVE on the CYBER 205 is given in 


Figure 7. The file containing this JCL can be permanently saved on the user's 


file space by typing the following commands: 
GET , RCP205J/UN=CER@Q2 


SAVE, RCP205J 
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This JCL for running RCPWAVE on the CYBER 205 must first pass through the 
CYBER 865 front-end computer. The first section of lines in the JCL deals 
with preparation of the program RCPWAVE on the front end and the subsequent 
relay to the CYBER 205. All JCL is described in detail in Appendix E. 


Submitting the CYBER 205 Batch Job 


73. The procedure for submitting a batch job for processing on the 
CYBER 205 is identical to that used for the CYBER 865. The following command 
sends the JCL shown in Figure 7 to the front-end CYBER 865 and ultimately to 
the CYBER 205: 


SUBMIT, RCP205J,T 
Computer response to the SUBMIT command informs the user of the computer gen- 
erated job name (AKQBRUD in this case). 


74. Status of the job sent to the front-end CYBER 865 is checked in the 


same manner as was presented before, using the ENQUIRE command: 
ENQUIRE, JN=RUD 

When the job AKQBRUD is in the PRINT-TERMINAL QUEUE, as shown below, 
AKQBRUD IN PRINT - TERMINAL QUEUE 


then a job should have been routed to the CYBER 205. Whether or not this oc- 


curred can be determined by checking the diary of events by typing 
OLD, RCPDAYF 
EST: 
or by typing 
QGET , RUD 
XEDIT , AKQBRUD 
and then using the editor to locate the diary of events in this file. The 


dayfile should contain the following lines, including the seven-letter jobname 
given to the file SENDJOB by the computer: 


13.40.08.SUBMIT,SENDJOB,T. 


13.40.08. SUBMIT COMPLETE. JOBNAME IS AKQBRUZ. 
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If this line is not found, errors have occurred in the batch job processing; 
and the job has not been relayed to the CYBER 205 machine. 

75. Assuming the job has been successfully sent from the front end to 
the CYBER 205, the status of the CYBER 205 job can be checked with the follow- 


ing command: 
LINK, ENQUIRE, JN=JOB205 
The computer will respond with 
STATUS FOR CER@Q2/KOE 


DATE TIME USER SYSTEM FILE FUE 
MMDD HHMM JOB NAME JOB NAME TYPE STATUS 


0326 1547 JOB205 AKQBRUZ TO ROUTING INITIATED TO ADY 
0326 1547 JOB205 JOB2RVB TO ARRIVED AT ADY 
LINK COMPLETE. 


When the job status reads 


STATUS FOR CER@Q2/KOE 


DATE TIME USER SYSTEM ETERS my Pale: 

MMDD HHMM JOB NAME JOB NAME TYPE STATUS 

0326 1547 JOB205 AKQBRUZ TO ROUTING INITIATED TO ADY 
0326 1547 JOB205 JOB2RVB TO ARRIVED AT ADY 

0326 1555 JOB205 AKQBRUZ WT OUTPUT AVAILABLE AT KOE/T 


then the batch job has been completed, and output is available to the user. 
76. In order to get a compiled version of RCPWAVE and a dayfile con- 
taining execution information and a summary of actual cost information from 


the CYBER 205, the following QGET command must be used: 
QGET , RUZ 
XEDIT , AKQBRUZ 


The text editor can then be used to locate and print the desired information. 
Notice that the output file RCPOTPT is absent from the CYBER 205 JCL. Cur- 
rently it is not possible to have this information saved as a permanent file 
within the JCL as was done on the CYBER 865. 

77. All other permanent files used or created by the JCL can be ac- 
cessed and viewed in any of the three ways given previously. For example, 


file RCPPRNT could be viewed by typing 


Wy 


GET, RCPPRNT 
XEDIT, RCPPRNT 
or 
OLD, RCPPRNT 
EIST 
or 
GET, RCPPRNT 


ROUTE, RCPPRNT , DC=PR, UN=SC 


where SC is the site code of the device which is to receive the printed file. 


Input Files 


78. The model source code RCPWAVE is the first file retrieved in the 
course of JCL processing. Both JCL files RCP865J and RCP205J obtain the 
source code from the permanent file space assigned to the user number CER@Q2. 
Users are discouraged from copying the source code onto their file space, mak- 
ing their own permanent modifications, and running their own model version. 
CERC personnel will support only the source code residing on this user number. 

79. A compiler-generated listing of the entire source code, which in- 
cludes a main program and nine subroutines, is provided in Appendix F. Users 
are encouraged to read comments within the program listing which explain im- 
portant variables and input parameters as well as those comments which de- 
scribe what is being performed in each subroutine. Another important feature 
contained in the listing is the line identifier to the far right of each 


FORTRAN statement. The following example of the sixth statement, 


PROGRAM RCPWAVE (INPUT,OUTPUT, TAPE7, TAPE6, TAPE8, TAPE3) MAIN 2 


has as its identifier "MAIN 2" . "MAIN" denotes the main program, and "2" 
denotes the second line. The identifier associated with each subroutine is 
the subroutine name itself, and, again, the numbers correspond to line numbers 
within each subroutine. These identifiers are required when making changes to 
the model source code using the "Update" method and ensure correct placement 


of changes. 
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80. Three additional input files are accessed by RCP865J and RCP205J. 
They must reside on the user's permanent file space. These files contain in- 
formation describing grid characteristics and wave conditions to be simulated, 
data controlling the amount of printed output, bottom bathymetry for each grid 
cell, and any changes to the model source code. The general format of each of 
these files is independent of the CDC machine being used to execute the wave 
model. However, one additional source code modification must be made when 
running on the CYBER 205. It will be discussed later in this section. In the 
following text, all three input files will be addressed using the names given 
to them in Figures 5 and 7--RCPDATA, RCPDEPT, and RCPUPDT. 

81. RCPDATA contains grid information, specifically the size and number 
of cells in the grid. It contains deepwater wave characteristics (height, pe- 
riod, and direction) describing wave conditions to be simulated, and it con- 
tains data which control the volume of printed output. These, along with all 
other input parameters, are described below. Following each parameter, in 
parentheses, is its corresponding value taken from the sample input data file 
shown in Figure 8. 

FIELD RESEARCH FACILITY (DUCK , NC) EXAMPLE 
rh) 12.0 50 24.0 9,800 1.0 2 0.90 
2.00 12.00 20.00 


1.50 8.00 -35.00 
dy Ue eet 4) 


Figure 8. Sample input data file for a 
wave model application 


82. A detailed description of the input parameters follows: 
a. CARD 1 - FORMAT (7A10). 


LTITLE - A 70-character array describing the model application 
being considered. 


(FIELD RESEARCH FACILITY (DUCK, NC) EXAMPLE) 


b. CARD 2 - FORMAT (15, F10.4, 15, OL AMO Sen Mn Oe) ic 
M - Number of grid cells in the x-direction (75) 
DX - Cell size in the x-direction (12.0) 
N - Number of cells in the y-direction (50) 
DY - Cell size in the y-direction (24.0) 
G - Gravitational acceleration constant (9.800) 


Note: The parameter G controls the units of all input 
and output variables. The choice shown above (9.800) 
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|Q. 
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requires that all input be in the kilograms, metres, 

and seconds unit system. All output will also be in 

these units. Other options for G are 32.2 (pounds, 
feet, and seconds) and 980.0 (grams, centimetres, and 
seconds). 


CNTRANG - Angle between bathymetric contours and the (0.0) 
y axis 


Note: If the bottom contours deviate from a line parallel to 
the y axis by a well defined angle, inclusion of this 
variable may result in faster iterative convergence to- 
ward a solution. The angle convention adopted for 
CNTRANG is shown in Figure 2. All angles are measured 
in degrees. If a well defined angle does not exist, or 
if there is doubt, set this parameter to zero. 


NCASES - Number of individual wave conditions to be (2) 
simulated 

DLEVEL - Constant added to the entire water depth (0.90) 
matrix which adjusts the water level datum 

CARDS 3a, 3b, ... - FORMAT (3F10.3). 

There must be one card for each wave condition considered. 

HDEEP - Deepwater wave height (2.00 for case 1, 1.50 for 
case 2) 

TDEEP - Wave period (12.00 for case 1, 8.00 for 
case 2) 

ZDEEP - Deepwater wave angle (20.00 for case 1, -35.00 for 
case 2) 


Note: See Figure 2 for the wave angle convention. All wave 
angles in both input and output files are given in 


degrees. 

CARD 4 - FORMAT (515). 

IP1 - Wave model output starts at this "I" row number (1) 

IP2 - Wave model output ends at this "I" row number (50) 

INC - Row by row output is incremented by this number (1) 
of lines 

JP1 - Wave model output starts at this "J" column Gan) 
number 

JP2 - Wave model output ends at this "J" column (41) 
number 


Note: Column by column output is always incremented by one. 
CARD 5 - FORMAT (215). 


IREF - Flag for instructing the model to consider (1) 
refraction only or combined refraction and 
diffraction 


= 0 refraction only 
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= 1 combined refraction and diffraction 
Note: Users should set this parameter to one. 


IGRID - Flag for indicating what type of grid system is 
being used (omitted) 


Note: This must be set to zero or omitted. 


83. The model can be run using either constant, or variably sized, 
rectangular grid cells. This model documentation discusses only applications 
involving grids with constant sized cells. As stated earlier, technology for 
creating variable sized grids which are compatible with RCPWAVE does exist at 
CERC. 

84. Bottom bathymetry for each grid cell must be supplied to the model. 
The source code is written such that bathymetric data are read from the file 
RCPDEPT in the following format. Depths along the entire row I=1 (from J=1 
to J=N) are read using a 10F8.2 format. Depths along row I=2 are read next, 
and the procedure is repeated until depths along the offshore row I=M have 
been read. An example of a bathymetric grid with 10 cells in the x-direction 
(I coordinate) and 24 cells in the y-direction (J coordinate) is shown in 
Figure 9. A depth file containing this data set, and written in 10F8.2 for- 
mat, is also shown. It is a trivial matter to change this section of the 
source code which reads the bathymetry (in subroutine DEPTH) to some other, 
more convenient format. Such a modification will be illustrated later in this 
section. 

85. The wave model considers only positive, nonzero water depths in its 
computations. All grid cells which occupy dry land can be designated by nega- 
tive depths or zeros in the input file, but the model will internally assign a 
small water depth of 1.0 ft (or a metric equivalent) to each of these cells. 
Depths in water cells that are less than 1.0 ft will be increased to this min- 
imum value. Bathymetric input data are modified within the program in one 
other way. Along each row, depths in columns J=1 and J=N are replaced by the 
depths in columns J=2 and J=N-1, respectively. These changes are consistent 
with the lateral boundary conditions used in the model. 

86. The third input file obtained by the JCL is the file RCPUPDT which 
contains necessary or desired corrections to the source code via the "Update" 
method. Figure 10 shows an example of an "Update" file. It is used to illus- 
trate certain points pertaining to this correction method. The first state- 


ment in the file 
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*ID MODS 


identifies all subsequent changes to be made with the identifier "MODS." 


Figure 9. 
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Sample bathymetry grid and input file for a wave model application 


line must be the first in any update file like RCPUPDT. The next two lines, 


*D PARAM. 3 
and 
PARAMETER (1Q=75,JQ=50) 


delete the line identified by "PARAM.3" everywhere in the program and replace 
it with the new PARAMETER statement. This action changes the limits on all 

arrays within the program to match the number of cells in the grid mesh. The 
quantity IQ must be set greater than or equal to the number of rows, M , 

and JQ must be greater than or equal to the number of columns, N . These 
first three "Update" lines are the only ones required for model applications 
run on the CYBER 865 machine. For jobs sent to the CYBER 205, two additional 


lines must be included, 


*D MAIN.2 
PROGRAM RCPWAVE (INPUT, OUTPUT, TAPE7=TAPE7 , TAPE6=TAPE6 , TAPE8=TAPE8 


* TAPE3=TAPE3) 


These are required because the PROGRAM statement format is different for each 


machine's compiler. They are shown also in Figure 10. 


ID MODS 
D PARAN. 3 
PARAMETER (1Q=75, JQ=50) 
D DEPTH. 23, DEPTH. 26 
DO 14 J=1,N 
READ (8,60) (D(1, J) > T=1.M) 
DO 315 I=1,M 
DT, J)=-DB(1, J) 
315 CONTINUE 
60 FORMAT (SX/ (15F8. 2) ) 
14 CONTINUE 
D MAIN. 203 
DNULT=10.0 
D MAIN.2 
PROGRAM RCPWAVE (INPUT OUTPUT TAPE7=TAPE7, TAPES=TAPEG, TAPEG=TAPES 
*, TAPES=TAPE3) 


Figure 10. Sample "Update" file for 
a wave model application 


87. The "Update" modification, 
*D DEPTH.23,DEPTH.26 
DO 14 J=1,N 


READ(8,60) (D(I,J),I=M) 
DO 315 I=1,M 
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Dera “DCUea) 
315 CONTINUE 
60 FORMAT(5X/(15F8.2)) 
14 CONTINUE 


and the seven subsequent lines of FORTRAN code delete the statements identi- 
fied by "DEPTH.23" through "DEPTH.26" in subroutine DEPTH and replace them 
with the seven new lines. These modifications alter the format used to read 
bottom bathymetry data; and, since the water depths on the file happen to be 
negative numbers, a Sign change is programmed into the procedure. The 
"Update" correction 


*D MAIN.203 
DMULT=10.0 


replaces the line in the main program identified by "MAIN.203" with the sub- 
sequent new line of code. This correction changes the accuracy of the depth 
matrix in the printed output from whole metres (or other units) to tenths of 
metres (or other units). As it is written, the model prints out water depths 
to the nearest whole foot (metre or centimetre), wave heights to the nearest 
tenth foot (metre or centimetre), and wave angles to the nearest degree. 
Changes in output accuracy must be made via "Update" corrections. Users are 
referred to the statements "MAIN.194" through "MAIN.206" in Appendix F for 
further explanation. 

88. These "Update" corrections have illustrated the use of the DELETE 
(*D) command which deletes line(s) of text. Other frequently used "Update" 
commands are the INSERT (*I) commmand, which inserts text after a certain 
line, and the BEFORE (*B) command which inserts text before a line. Users 
interested in making source code changes should first consult the "Update" 


manual referenced earlier. 


Output Files 


89. Two permanent output files are generated during the course of model 
execution on the CYBER 205. These same two files plus one additional file are 
generated during the course of execution on the front end (CYBER 865). Output 
file names given in Figure 7 will be used in this discussion. Information 


contained on the printed output file RCPPRNT is independent of the choice of 


Dil: 


machines, but the diary of events RCPDAYF is machine dependent. The file 
RCPOTPT is generated only by the front-end machine. At some point in the 
future, RCPOTPT will be generated also by the CYBER 205. 

90. The RCPPRNT file summarizes all model input such as grid charac- 
teristics, wave conditions considered, and bathymetry for that window of 
the grid mesh defined by the parameters IP1, IP2, INC, JP1, and JP2. For 
each wave condition simulated, a listing of the following wave characteristics 
is printed for the two-dimensional field defined by the above parameters: 

(a) breaker indices defining the location of the surf zone (zeros denote 

cells where no breaking is occurring), (b) wave angles, (c) wave heights, and 
(d) the magnitude of the wave phase function gradient (in the absence of ap- 
preciable diffractive effects this is approximately the wave number, e.g., it 
is a function of the wave length). Specific examples of this file will not be 
presented here but rather in Part V, the sample application section. 

91. For jobs run on the CYBER 865, a short file (RCPDAYF) is created 
which summarizes the steps performed during the course of JCL processing. Any 
errors encountered during compilation and/or execution are indicated on this 
file. Estimated costs are given also, but these are not the costs which are 
actually billed to the user's account. Even though this "dayfile" indicates 
any errors which have caused abnormal job termination, it provides little 
detail for locating the problem. These details are provided on the file 
RCPOTPT. 

92. This rather large file contains a compiled listing of the model, 
any compilation errors which may have occurred, and a short variable map. Be- 
cause of its length, an example is not printed. Any of these output files can 
be viewed using one of the three methods mentioned previously. It should be 
noted that the job output returned to the users terminal with the T parameter 
on the SUBMIT command is very useful and can be viewed using the QGET command. 
This file contains: (a) the entire RCPDAYF file, (b) actual costs (including 
the Corps of Engineers discount), and (c) the entire RCPOTPT file. 

93. For CYBER 205 model execution the "dayfile," RCPDAYF, contains only 
those steps actually processed on the front-end CYBER 865 machine. The QGET 
command must be used to obtain job output from the CYBER 205. This output in- 
cludes: (a) a file equivalent to the front-end RCPOTPT file (compiled listing 
and any compilation errors), (b) the diary of events which occurred during 


CYBER 205 processing, and (c) computer usage and actual cost information. 
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PART V: MODEL APPLICATIONS 


General Comments 


94. General guidance concerning RCPWAVE applications is provided in 
this section. Two model applications are given to demonstrate the different 
scales of prototype problems for which the model can be applied successfully. 
The next two sections explain each example in greater detail. Listings of JCL 
files used to run each simulation are provided, as are listings of input files 
and printed output. Comparisons between model execution costs incurred on 
both CDC computers are given in the last section. 

95. Whenever RCPWAVE is a candidate for solving a wave propagation 
problem, the area of interest is well defined, and wave conditions to be simu- 
lated are known. The user must create a finite difference grid system encom- 
passing the area. Grid characteristics determine the success of the applica- 
tion, so they should be given considerable thought. Two factors completely 
define the grid, orientation of the axes, and cell resolution. 

96. Users should construct the grid system so that the y-axis runs as 
parallel to the coastline as possible. This causes the x-axis to be directed 
offshore and probably somewhat perpendicular to the bottom contours. The lat- 
eral boundary conditions used in the model are most accurate when this axis 
configuration is adopted. The boundary conditions assume that the variation 
of the bathymetry, and wave parameters, in the y-direction is small. The as- 
sumption is most accurate when the y-direction is essentially the longshore 
direction. Ideally this axis configuration also results in the offshore ba- 
thymetry having fairly straight and parallel contours and the offshore rows of 
cells having the greatest depths. The procedure of applying Snell's law to 
propagate wave information from deep water to the seaward boundary of the 
model is more accurate if these conditions exist. 

97. Waves will generally attempt to advance in directions parallel to 
the x-axis as a result of the recommended orientation. The model tends to 
produce erroneous results if computed waves attempt to propagate parallel to 
the y-axis (wave angles near plus or minus 90 deg). Even with this axis con- 
figuration, problems may arise due to very irregular bottom bathymetry and/or 
very oblique wave incidence. These errors manifest themselves via large local 


Wave angles and occurrences of very large and very small wave heights in 
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alternating cells. An artificial limit on the absolute value of the wave an- 
gle, 86 deg, is written into the model. Problems of this nature can sometimes 
-be eliminated by using finer grid resolution. 

98. Cell size is the second factor defining the grid system. Accuracy 
of model solutions computed using finite difference approximations is depen- 
dent upon cell size; smaller cells result in less error. Diffractive effects 
can become important on spatial scales with orders of magnitude smaller than 
the wave lengths being considered. In modeled regions of interest where dif- 
fraction is important (very complicated bathymetry) cell sizes may need to be 
some fraction of the wavelength in order to simulate wave propagation accu- 
rately. Fine resolution is not a requirement for model applications where the 
bathymetry is very regular. The physical processes of importance, whether it 
be refraction or a combination of refraction and diffraction, dictate cell 
sizes which should be used. 

99. The examples in the next two sections represent two different kinds 
of problems. Example I deals with wave propagation over a rather small region 
with complex bathymetry. Diffraction was presumed to be an important physical 
process; therefore, finer grid resolution was used in an attempt to represent 
this effect more accurately. Example II illustrates the applicability of 
RCPWAVE for propagating waves over a very large area dominated by very regular 
bathymetry. Assuming that refractive effects would dominate wave transforma- 
tion throughout the region, cell sizes were increased to a value on the same 
order of magnitude as the wave length. These two examples show the diversifi- 


cation of problems which can be addressed using the model. 


Example I: FRF Pier, Duck, North Carolina 


100. Bathymetry in the vicinity of CERC's research pier at Duck, North 
Carolina, is quite complex as seen in Figures 11 and B2. Figure 11 was 
plotted using similar bathymetric data digitized onto the rectangular grid 
system comprised of 75 cells in the offshore (x) direction and 50 cells in 
the longshore (y) direction. Again, cell sizes are 12 and 24 m in the x- and 
y-directions, respectively. All depths on the bathymetric file are given in 
metres below MSL. A constant tidal elevation of 0.90 m is added to the 
bathymetric data. Two wave conditions are considered: (a) 2-m-deep water 


wave height, 12-sec wave period, and +20 deg deepwater incident angle (refer 
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Y-AXIS 


MODEL OUTPUT 
FOR THIS REGION 


DUCK PIER, NORTH CAROLINA BATHYMETRY 


X- AXIS DEPTHS ARE METERS BELOW MSL 


Figure 11. Bathymetry used in the FRF pier, Duck, North Carolina 
sample application 


to the angle convention in Figure 2) and (b) 1.5-m deepwater wave height, 
8-sec period, and -35 deg deepwater incident angle. 

101. All files associated with this example are listed in Appendix G. 
Files are arranged in the following sequence: 


a. JCL file to run this application on the CYBER 865, DUCK865 
(Figure G1) 


b. JCL file to run this application on the CYBER 205, DUCK205 
(Figure G2) 


ec. Input data file , DUCKDAT 
(Figure G3) 


d. Source code correction file for the CYBER 865 run, DUCKUPD 
(Figure G4) 
e. Source code correction file for the CYBER 205 run, DUCKUP2 


(Figure G5) 
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f. A sample of the bathymetric data file , DUCKDEP 
(Figure G6) 
g. Printed output file , RCPPRNT 


(Figure G7) 


Two-dimensional fields of total water depth, breaker index, wave angle, wave 
height, and wave number (more precisely, the wave phase gradient) are printed 
for a portion of the grid extending from rows I=1 to I=50 and from columns 
J=11 to J=41. Figure 11 shows the extent of this region relative to the en- 
tire grid system. 


Example II: Homer Spit, Alaska 


102. Bathymetry for this application (see Figure 12) is digitized 
onto a grid with 96 rows in the offshore (x) direction and 83 cells in the 
longshore (y) direction. Cell dimensions are approximately 417 ft in the 
x-direction and 833 ft in the y-direction (about 10 times the resolution used 


in the Duck, North Carolina, examples. One wave condition is considered in 
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Figure 12. Bathymetry used in the Homer Spit, Alaska, sample application 
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the simulation. The deepwater wave height is 8 ft, the period is 10 sec, and 
the deepwater wave angle is -45 deg (see Figure 2 for the wave angle conven- 
tion). Bathymetry is given in fathoms below mean low water so depths are con- 
verted to feet via a correction to the source code. No water level datum ad- 
justment is applied. 

103. All files associated with this example are listed in Appendix H. 
Files are arranged in the following sequence: 


a. JCL file to run this application on the CYBER 865, HOMR865 
(Figure H1) 


b. JCL file to run this application on the CYBER 205, HOMR205 
(Figure H2) 


ec. Input data file , HOMRDAT 
(Figure H3) 


d. Source code correction file for the CYBER 865 run, HOMRUPD 
(Figure H4) 


e. Source code correction file for the CYBER 205 run, HOMRUP2 
(Figure H5) 

f. A sample of the bathymetric data file , HOMRDEP 
(Figure H6) 

g. Printed output file , RCPPRNT 


(Figure H7) 


Two-dimensional fields of water depth and wave parameters are given also in 


Appendix H for the subgrid region denoted in Figure 12. 


Cost Comparisons 


104. The following tabulation shows cost comparisons between simula- 
tions run on both CDC machines. Costs in dollars are shown for Examples I and 


II which represent the FRF pier and Homer Spit, respectively. 


CYBER 865 CYBER 205 
Priority Cost _ in Dollars Cost in Dollars 
Example I Py 2.24 6.09 
P3 0.86 4.74 
P2 Oni 1.68 
Example II PY 2051 6.40 
P3 0.96 4.98 
P2 0.19 Bed 


Dilh 


These are actual costs to the user and include the Corps' fiscal year 1985 
(FY 85) discount of 45 percent. Example I was run on a 75 by 50 grid. Two 
wave conditions were simulated. One wave condition was considered in Exam- 
ple II, but the grid was larger, 96 by 83. 

105. The CYBER 865 is a cheaper, easier computer to work on. It does 
have central memory limitations. Comparisons show that costs incurred during 
simulations using the CYBER 205 are not excessive. Users should not sacrifice 


modeling accuracy simply to avoid using the CYBER 205. 
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PART VI: GRID INTERPOLATION PROGRAM (INTPRCP) 
General Comments 


106. Bathymetry at every cell of the finite difference grid mesh is re- 
quired as input into RCPWAVE. Since digitization of bathymetric data from 
survey or nautical charts is both tedious and time consuming, it is desirable 
to limit the number of times one must perform this task. The grid interpola- 
tion program INTPRCP provides a method for determining bathymetry on a "new" 
grid using depth data from an "old" grid. The origin of the new grid can be 
rotated and/or translated relative to the old grid origin. This capability 
for interpolation ensures that the bulk of the digitization is done only once. 
The program is useful when: 


a. A coarse grid containing bathymetry is available but finer 
resolution is desired or required. 


b. The axis orientation of the available grid is different from 
the desired orientation. 


c. The region to be modeled is a subset of the available grid. 


A compiled listing of the program is given in Appendix I. 

107. The program uses a fairly simple interpolational scheme. For each 
new grid cell, a search is done for the nearest old depth value in each of the 
four 90 deg quadrants relative to the new cell center. A weighted average of 
these four values is used to compute a new depth. The weighting functions are 
based on relative distances from the new grid center to the cell centers de- 
fining the positions of the four old depth values. 

108. JCL is available to run the program on the CDC CYBER 205 machine 
only. The next section briefly describes the procedures performed in the 
course of job execution. It also explains how an interpolational job is sub- 
mitted. The third section documents, in detail, all required input files and 
includes some examples. It also describes the output files created during 
program execution. A complete example illustrating the use of INTPRCP is 


presented in the last section. 


Executing INTPRCP on the CDC Computing System 


109. The interpolation program is run on the CYBER 205 machine in a 
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remote batch mode just like the wave model. This machine was chosen instead 
of the CYBER 865 because it has more available central memory. As a result, 
users are not limited by the number of grid cells contained in either the 
old or new grids. A CYBER 205 user number and password must be obtained be- 
fore attempting to run INTPRCP. The direct access file AF205, mentioned in 
Part IV, must exist on the user's file space. This file contains validation 
information for the user's CYBER 205 and CYBER 865 accounts. Refer to Part IV 
for details concerning these requirements. 

110. The JCL file called INTRCPJ, used to run INTPRCP, is shown in Fig- 
ure 13. A copy of this file can be obtained by logging into the CDC computing 
system and typing 


GET , INTRCPJ/UN=CEROQ2 


This action creates a local file called INTRCPJ on the user's work space. To 


save the file permanently, type the following command: 


SAVE, INTRCPJ 


The JCL file for executing the interpolation program on the CYBER 205 must 
first pass through the front-end machine. The first section of JCL commands 
deals with compilation of the program on the front end and the subsequent re- 
lay to the CYBER 205. Commands in the CYBER 205 JCL instruct the computer to 
gather old and new grid information, retrieve bathymetry data for the old 
grid, compile and execute the program, and save the new grid bathymetry. The 
"COMMENT" section of the JCL briefly describes required input data files and 
output files generated in the course of program execution. Before submitting 
this JCL file, the user must make a few changes. Users must replace the "<205 
USERNAME>," "<205 PASSWORD>," and "<CYBER 865 USERNAME>" with the passwords 
and user numbers associated with their account. 

111. The front-end machine is given the entire batch of commands con- 


tained in INTRCPJ by using the SUBMIT command in the following way: 

SUBMIT, INTRCPJ,T 
The T parameter instructs the computer to return certain output to the user's 
terminal. The procedure for checking job status and accessing output is iden- 


tical to that presented in Part IV. No further discussion addressing these 


points is given. 
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Input and Output Files 


112. The interpolation program reads input data from two files, INTDEPO 
and INTPGRD. Old grid bathymetry is contained in INTDEPO. The procedure for 
reading bathymetry into this program is the same as that used in RCPWAVE. The 
row of depths corresponding to I=1 are read first in "10F8.2" format. A row 
includes values from J=1 to J=N. The row-by-row procedure is repeated until 
the offshore row (I=M) has been read. An example of a depth file in this for- 
mat was shown in Figure 9. 

113. Input data contained in INTPGRD are geometric parameters which de- 
fine the grid characteristics of the old and new grids. The format used by 
the program to read this information and a description of each parameter are 
given here. The value for each input parameter, taken from the sample input 


file shown in Figure 14, is also given in parentheses after the description. 


35 30) «650 40. 10.0 0.5 0.9 
Ono 

0.20 0.20 

0.14 0.15 


Figure 14. Sample input data file 
for grid interpolation 


a. CARD 1 - FORMAT (415, 3F10.5). 
M - Number of grid cells in the x-direction for (35) 
the old grid 
N - Number of grid cells in the y-direction for (30) 
the old grid 
N2 - Number of grid cells in the x-direction for (50) 


the new grid 


M2 - Number of grid cells in the y-direction for (40) 
the new grid 


ANG - Angle of rotation measured from the old grid (10.0) 
to the new grid (positive angles are meas- 
ured counterclockwise, negative angles are 
measured clockwise) 


XSHIFT - Offset from the old grid to the new grid in (0.5) 
the x-direction (positive offsets indicate 
that the new grid origin is below the old 
origin, negative offsets indicate that the 
new grid origin is above the old grid origin) 
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XSHIFT 


In 


10 


YSHIFT - Offset from the old grid to the new grid in 
the y-direction (positive offsets indicate 
that the new grid origin is to the right of 
the old origin, negative offsets indicate 
that the new grid origin is to the left of the 
old grid origin) 

Note: The conventions for ANG, XSHIFT, and YSHIFT are 

shown in Figure 15. 


YSHIFT OLD GRID 
SYSTEM 


1 
! 
! 
| 
| 
1 NEW GRID SYSTEM 
1 


Xoo 


XNEW 


Figure 15. "Old" and "new" grid conventions used 
in the interpolation program 


CARD 2 - FORMAT (215). 

GRDTYP1 - Old grid type (must be set to zero) 
GRDTYP2 - New grid type (must be set to zero) 
CARD 3 - FORMAT (2F10.5). 

DX - Old grid cell size in x-direction 


DY - Old grid cell size in the y-direction 
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(0.5) 


(0) 
(0) 


(0.2) 
(0.2) 


d. CARD 4 - FORMAT (2F10.5). 
DX2 - New grid cell size in the x-direction (0.14) 
DY2 - New grid cell size in the y-direction (0.15) 


Note: The parameters DX, DY, DX2, and DY2 must all 
have consistent units. 

114. In addition to bathymetry and grid information, one other input 
file is required. This is the file INTPUPD which contains any source code 
changes and must include information specifying the number of cells contained 
in both the old and new grids. An example of this file is shown in Figure 16. 
Variables in the PARAMETER statements are defined here; and their values, 


taken from the example, are shown in parentheses: 


IQ - Number of cells in the x-direction for the old grid (35) 
JQ - Number of cells in the y-direction for the old grid (30) 
IQ2 - Number of cells in the x-direction for the new grid (50) 
JQ2 - Number of cells in the y-direction for the new grid (40) 
KQ - Larger of the two products IQ*JQ or 1Q2*JQ2 (2000) 


Modifications to change the procedure or format for reading the bathymetry 
into the program should be included in this file. Figure 16 shows an example 
of such changes. Here, a format different from the one in the program was 
desired. The section of source code which reads old grid bathymetry is con- 
tained in the lines designated by "INTERPO.15" through "INTERPO.18" in Appen- 
Gixs Ii: 


41D MODS 
4D PARAH.3, PARAM. 5 
PAPAMETER(1Q=35, JQ=30) 
PARAMETER( 102=50 ,JQ2=40) 
PARAMETER (KQ=2000) 
4D INTERPO.18 
11 FORMAT(15F5.0) 


Figure 16. Sample 

"Update" file for 

grid interpolation 

115. There are three output files created as a result of running the 

program. The file INTPDAY is a dayfile containing a diary of events which 
occurred on the front-end CYBER 865 machine. This is saved as a permanent 
file on the user's file space. A file called INTPOUT is created on the 
CYBER 205 and is saved permanently. It contains a compile listing, execution 


information, and printed output created during the run. The printed output 
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includes a listing of both old and new grids and the grid characteristics of 
each. The file INTDEPN contains the bathymetry for the new grid. This file 
is written using the same procedure and format as those used for reading the 
old depth data (row by row in 10F8.2 format). 


Grid Interpolation Example 


116. The example problem involves interpolating bathymetry from an old 
grid with square cells to a new grid with smaller, rectangular-shaped cells. 
The new grid system is translated and rotated relative to the old grid. Orig- 
inal grid cell dimensions are 0.2 and 0.2 in the x- and y-directions, respec- 
tively. The new grid requires 50 cells in the x-direction, each 0.14 in 
length, and 40 cells in the y-direction, each 0.15 in length. The new grid is 
shifted by 0.5 in both coordinate directions and rotated 10 deg in the coun- 
terclockwise direction relative to the old grid. The input files used in this 
example, INTPGRD, INTDEPO, and INTPUPD are shown in Appendix J. The order in 


which the files are listed is: 


a. Input data file , INTPGRD (Figure J1) 
b. Source code correction file , INTPUPD (Figure J2) 
c. Sample of the input bathymetric file  ,INTDEPO (Figure J3) 
d. Sample of the output bathymetric file ,INTDEPN (Figure J4) 


Figure 17 shows both grids and contour plots of the original and interpolated 


bathymetry. 
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y-NEW 


YOLD 


NEW 
GRID 


ORIGINAL BATHYMETRIC CONTOURS 
ET INTERPOLATED BATHYMETRIC CONTOURS 


Figure 17. Comparison between original and interpolated bathymetry 
for the sample grid interpolation application 
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PART VII: CONCLUSIONS AND RECOMMENDATIONS 


117. The model presented in this report is capable of simulating wave 
propagation over arbitrary, and potentially complex, bathymetry. The govern- 
ing equations are theoretically applicable to linear, monochromatic plane wave 
propagation only. Comparisons between model results and laboratory data 
showed good agreement, indicating that the model is quite accurate if the 
above conditions are satisfied. The results also showed the model's ability 
to simulate not only wave refraction but also bottom-induced diffraction. 
Comparisons between model results and field data indicated the model is capa- 
ble of accurately simulating the transformation of single-peaked, narrow- 
banded, wave spectra. Results from one test suggest potential model use for 
solving problems involving propagation of wider-banded spectra, if a clearly 
defined wave direction exists. This hypothesis has not been substantiated in 
this report. 

118. The wave breaking scheme incorporated into the model is a vast 
improvement over those wave propagation models which assume proportionality 
between wave height and water depth throughout the surf zone. Comparisons be- 
tween model results and laboratory data showed the model to be reasonably ac- 
curate for a variety of bottom profiles, including plane, stepped, and barred 
beaches. Since the breaking scheme assumes that all waves break at the same 
location, the model is most valid for monochromatic waves and very narrow- 
banded wave spectra. No field data were used to verify the wave breaking as- 
pects of the model. 

119. The user's manual portion of the documentation and the sample ap- 
plications presented show the model's ease of application. Users are reminded 
to pay close attention to comments made concerning relationships between model 
accuracy and grid resolution. The problems of interest and the perceived 
physical processes of importance should dictate cell sizes. Also, the user 
should remember assumptions inherent in the model's governing equations. 
Applications in which these assumptions are violated may yield erroneous 
results. 

120. The model provides a tool which can be used by field personnel to 
solve many types of wave problems in a quantitative way and at a very low 
cost. If all assumptions inherent in the model development are essentially 


met, very accurate results can be obtained. Even if the model is applied to 
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problems in which some of the model assumptions are violated, it may still 
yield informative results, at least qualitatively if not quantitatively. 

121. The JCL files used to execute the wave model and the interpolation 
program on the CDC computing system utilize commands and procedures which are 
currently operative (at the time the report is being written). CERC personnel 
are not responsible for any changes implemented by the CDC which may render 
these files inoperative. However, CERC personnel will attempt to maintain us- 


able JCL files with the same names as those given in the text. 


68 


REFERENCES 


Abbott, M. B. 1975. Computational Hydraulics, Elements of the Theory of Free 
Surface Flows, Fearon-Pitman, Belmont, Calif. 
Battjes, J. A., and Janssen, J. P. 1978. "Energy Loss and Setup Due to 
Breaking of Random Waves," Proceedings of the 16th International Conference on 
Coastal Engineering, American Society of Civil Engineers, pp 569-587. 
Berkhoff, J. C. W. 1972. "Computation of Combined Refraction-Diffraction," 
Proceedings of the 13th International Conference on Coastal Engineering, Amer- 
ican Society of Civil Engineers, Vol 1, pp 471-490. 

1976. "Mathematical Models for Simple Harmonic Linear Water 


Waves, Wave Diffraction and Refraction," Publication No. 1963, Delft Hydrau- 
lics Laboratory, Delft, The Netherlands. 


Berkhoff, J. C. W., Booy, N., and Radder, A. C. 1982. "Verification of Nu- 
merical Wave Propagation Models for Simple Harmonic Linear Water Waves," 


Coastal Engineering, Vol 6, pp 255-279. 
Booij, N. 1981. "Gravity Waves on Water with Non-Uniform Depth and Current," 


Report No. 81-1, Department of Civil Engineering, Delft University of Technol- 
ogy, Delft, The Netherlands. 


Candel, S. M. 1979. “Numerical Solution of Wave Scattering Problems in a 
Parabolic Approximation," Journal of Fluid Mechanics, Vol 90, Part 3, 
pp 467-507. 


Dally, W. R. 1980 (May). "A Numerical Model for Beach Profile Evaluation," 
Master's Thesis, University of Delaware, Newark, Del. 


Dally, W. R., Dean, R. G., and Dalrymple, R. A. 1984. "Modeling Wave Trans- 
formation in the Surf Zone," Miscellaneous Paper CERC-84-8, US Army Engineer 
Waterways Experiment Station, Vicksburg, Miss. 


Divoky, D., Le Méhauté, B., and Lin, A. 1970. "Breaking Waves on Gentle 


Slopes," Journal of Geophysical Research, Vol 75, No. 9, pp 1681-1692. 


Dobson, R. S. 1967. "Some Applications of a Digital Computer to Hydraulic 
Engineering Problems," Technical Report No. 80, Department of Civil Engineer- 
ing, Stanford University, Stanford, Calif. 


Dunham, J. W. 1951. "Refraction and Diffraction Diagrams," Proceedings of 
the 1st Conference on Coastal Engineering, Council on Wave Research, Engineer- 
ing Foundation, pp 33-49. 

Goda, Y. 1970. "A Synthesis of Breaker Indices," Transactions of the Japa- 
nese Society of Civil Engineers, Vol 2, Part 2. 

Harrison, W., and Wilson, W. S. 1964. "Development of a Method for Numerical 
Calculation of Wave Refraction," Technical Memorandum No. 6, US Army Engineer 
Waterways Experiment Station, Coastal Engineering Research Center, Vicksburg, 
Miss. 

Horikawa, K., and Kuo, C. T. 1966. "A Study of Wave Transformation Inside 


the Surf Zone," Proceedings of the 10th International Conference on Coastal 
Engineering, American Society of Civil Engineers, pp 217-233. 


69 


Houston, J. R. 1981. "Combined Refraction and Diffraction of Short Waves 
Using the Finite Element Method," Applied Ocean Research, Vol 3, No. 4, 
pp 163-170. 


Hubertz, J. M. 1981. "Prediction of Wave Refraction and Shoaling Using Two 
Numerical Models," Coastal Engineering Technical Aid No. 81-12, US Army En- 
gineer Waterways Experiment Station, Coastal Engineering Research Center, 
Vicksburg, Miss. 


1982. "Prediction of Nearshore Wave Transformation," Coastal En- 
gineering Technical Aid No. 82-7, US Army Engineer Waterways Experiment Sta- 
tion, Coastal Engineering Research Center, Vicksburg, Miss. 


Iwata, K., and Sawaragi, T. 1982. "Wave Deformation in the Surf Zone," Mem- 


oirs of the Faculty of Engineering, Vol 34, No. 2, Nagoya University, Nagoya, 
Japan, pp 239-283. 

Izumiya, T. 1984. "A Study of Wave and Wave-Induced Nearshore Currents in 
the Surf Zone," Ph. D. Dissertation, University of Tokyo, Tokyo, Japan. 


Johnson, J. W., O'Brien, M. P., and Issacs, J. D. 1948 (Jan). "Graphical 
Construction of Refraction Diagrams," HO No. 605, TR-2, US Naval Oceanographic 
Office, Washington, DC. 


Kirby, J. T. 1983. "Propagation of Weakly-Nonlinear Surface Gravity Waves in 
Regions of Varying Depth and Current," Research Report No. CE-83-37, Depart- 
ment of Civil Engineering, University of Delaware, Newark, Del. 


Komar, P. D., and Gaughan, M. K. 1972. "Airy Wave Theory and Breaker Height 
Prediction," Proceedings of the 13th International Conference on Coastal Engi- 
neering, American Society of Civil Engineers, pp 405-418. 

Le Méhauté, B., and Koh, R. C. Y. 1967. "On the Breaking of Waves Arriving 
at an Angle to the Shore," Journal of Hydraulic Research, Vol 5, No. 1, 

pp 67-88. 

Lozano, C., and Liu, P. L. F. 1980. "Refraction/Diffraction Model for Linear 
Surface Water Wave," Journal of Fluid Mechanics, Vol 101, Part 4, pp 705-720. 


McCowan, J. 1891. "On the Solitary Wave," Philosophical Magazine, 5th Se- 
ries, Vol 32, No 134, pp 45-58. 


Miche, R. 1944. "Mouvements Ondulatoires de la Mer en Profondeur Constante 
ou Decroissante," Annales des Ponts et Chaussees, Series 3, Issue 363, 
pp 25-78, 131-164, 270-292, and 369-406. 


Mizuguchi, M. 1980. "A Heuristic Model of Wave Height Distribution in the 


Surf Zone," Proceedings of the 17th International Conference on Coastal Engi- 
neering, American Society of Civil Engineers, pp 278-289. 


Noda, E. K., et al. 1974. "Nearshore Circulations Under Sea Breeze Condi- 
tions and Wave-Current Interactions in the Surf Zone," Technical Report No. 4, 
Tetra-Tech, Inc., Pasadena, Calif. 


Perlin, M., and Dean, R. G. 1983. "A Numerical Model to Simulate Sedi- 
ment Transport Within the Vicinity of Coastal Structures," Miscellaneous 
Paper 83-10, US Army Engineer Waterways Experiment Station, Coastal 
Engineering Research Center, Vicksburg, Miss. 


70 


Pierson, W. J., Neumann, G., and James, R. W. 1952. "Observing and Forecast- 
ing Ocean Waves," HO No. 603, US Naval Oceanographic Office, Washington, DC. 


Rabe, K. 1975. "The Delaware-Dobson Wave Refraction Model," Computer Pro- 
gramming Note No. 21, Environmental Prediction Research Facility, Naval Post- 
graduate School, Monterey, Calif. 


Radder, A. C. 1979. "On the Parabolic Equation Method for Water-Wave Propa- 
gation," Journal of Fluid Mechanics, Vol 95, Part I, pp 159-176. 


Sheng, Y. P., Segur, H., and Lewellen, W. S. 1978. "Application of a Spatial 
Smoothing Scheme to Control Short-Wave Numerical Oscillations," Technical Mem- 
orandum No. 78-8, Aeronautical Research Associates of Princeton, Princeton, 
New: 


Smith, R., and Sprinks, T. 1975. "Scattering of Surface Waves by a Conical 
Island," Journal of Fluid Mechanics, Vol 72, Part 2, pp 373-384. 


Thornton, E. B., and Guza, R. T. 1982. "Energy Saturation and Phase Speeds 


Measured on a Natural Beach," Journal of Geophysical Research, Vol 87, 

No. C12, pp 9499-9508. 

Tsay, T.-K., and Liu, P. L.-F. 1982. "Numerical Solution of Water-Wave Re- 
fraction and Diffraction Problems in the Parabolic Approximation," Journal of 
Geophysical Research, Vol 87, No. C10, pp 7932-7940. 


Weggel, J. R. 1972. "Maximum Breaker Height," Journal of the Waterways, Har- 
bors, and Coastal Engineering Division, Vol 78, No. WW4, pp 529-548. 

Whalin, R. W. 1971. "The Limit of Applicability of Linear Wave Refraction 
Theory in a Convergence Zone," Research Report H-71-3, US Army Engineer Water- 
ways Experiment Station, Coastal Engineering Research Center, Vicksburg, Miss. 


. 1972. "Wave Refraction Theory in a Convergence Zone," Proceed- 
ings of the 13th International Conference on Coastal Engineering, American So- 
ciety of Civil Engineers, Vol 1, pp 451-470. 

Williams, R. G., Darbyshire, J., and Holmes, P. 1980. "Wave Refraction and 
Diffraction in a Caustic Region: A Numerical Solution and Experimental Vali- 


dation," Proceedings of the Institute of Civil Engineering, Vol 69, Part 2, 
pp 635-649. 


a 


(aad te ia’ 
ce af sds 


APPENDIX A: VERIFICATION OF MODEL RESULTS USING THE 
ELLIPTICAL SHOAL LABORATORY TEST DATA 


WAVE DIRECTION 


a0 15 10 5 0 
Y COORDINATE, METERS 


Figure Al. Geometry used in the model simulation of the 
elliptical shoal case 
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Figure A2. Comparisons between model results and observed data 
for the elliptical shoal case (profiles 1, 2, and 3) 
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Figure A3. Comparisons between model results and observed data 
for the elliptical shoal case (profiles 4 and 5) 
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Figure A4. Comparisons between model results and observed 
data for the elliptical shoal case (profiles 6, 7, and 8) 
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APPENDIX B: VERIFICATION OF MODEL RESULTS USING FIELD 
RESEARCH FACILITY (FRF) PROTOTYPE DATA 
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Figure B1. 


Observed offshore wave spectra during field 


verification, Cases 1 through 6 


B2 


Y-AXIS 


OCT 16, 1982 BATHYMETRY 
CONTOURS ARE GIVEN IN METERS 
BELOW MEAN SEA LEVEL 


DISTORTED SCALE 


xX - AXIS 


Y-AXIS 


PIER 


OCT 27, 1982 BATHYMETRY 
CONTOURS ARE GIVEN IN METERS 
BELOW MEAN SEA LEVEL 


ee ah ln UE 


DISTORTED SCALE 


nn 


= 


X- AXIS ASS 


Figure B2. Bathymetry around the FRF pier measured 
on 16 October 1982 and 27 October 1982 
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Figure B3. Comparisons between model results and observed 
data for field verification, Cases 1, 2, and 3 
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Figure B4. Comparisons between model results 
and observed data for field verification, 
Cases 4, 5, and 6 
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APPENDIX C: VERIFICATION OF MODEL RESULTS USING BREAKING 
WAVE DATA COLLECTED IN LABORATORY EXPERIMENTS 
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Figure C1. Comparisons between model results and 
observed data for wave breaking verification, 
Test 1 (surf zone transformation only) 
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Figure C2. Comparisons between model results and ob- 
served data for wave breaking verification, Test 2 
(surf zone transformation only) 
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Figure C3. Comparisons between model results and ob- 
served data for wave breaking verification, Test 3 
(surf zone transformation only) 
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Figure C4. Comparisons between model results and ob- 
served data for wave breaking verification, Test 4 
(surf zone transformation only) 


C3 


WAVE HEIGHT (CM) 


DEPTH (CM) 


WAVE HEIGHT (CM) 


DEPTH (CM) 


@ 420 820 1200 1680 2020 2400 
DISTANCE OFFSHORE (CM) 
LEGEND 
= REP WAVER MODEL 
2... EXPERIMENTAL DATA HORIKAWA LAB TEST 5 
iS BOON PRO hole e (BREAK POINT MATCH) 


Figure C5. Comparisons between model results and ob- 
served data for wave breaking verification, Test 5 
(surf zone transformation only) 
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Figure C6. Comparisons between model results and ob- 
served data for wave breaking verification, Test 6 
(surf zone transformation only) 
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Figure C7. Comparisons between model results and ob- 
served data for wave breaking verification, Test 7 
(surf zone transformation only) 
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Figure C8. Comparisons between model results and ob- 
served data for wave breaking verification, Test 8 
(surf zone transformation only) 
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Figure C9. Comparisons between model results and ob- 
served data for wave breaking verification, Test 9 
(surf zone transformation only) 
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Figure C10. Comparisons between model results and ob- 


served data for wave breaking verification, Test 5 
(incipient breaking and surf zone transformation) 
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Figure C11. Comparisons between model results and ob- 
served data for wave breaking verification, Test 7 
(incipient breaking and surf zone transformation) 


WAVE HEIGHT (CM) 


DEPTH (CM) 


+ - + ae + | 


* 
1080 288 300 400 580 600 788 
DISTANCE OFFSHORE (CM) 


LEGEND 
oS RERWAVE MODEL 
.+.-. EXPERIMENTAL DATA LZUM EVA EABS TES 6 
oS — SO mIOMPROR TIE 


Figure C12. Comparisons between model results and ob- 
served data for wave breaking verification, Test 8 
(incipient breaking and surf zone transformation) 
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Figure C13. Comparisons between model results and ob- 
served data for wave breaking verification, Test 9 
(incipient breaking and surf zone transformation) 
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APPENDIX D: LINE-BY-LINE DESCRIPTION OF JOB CONTROL LANGUAGE 
FOR EXECUTING RCPWAVE ON THE CYBER 865 COMPUTER 


(1) 


(2) 


(3) 


(4) 


/ JOB 


JOB , P3,T1000,CM37000. 


/USER 


/ CHARGE 


**Tmportant** 


The job card tells the computer that this 
is a batch job. 


The resource card informs the computer of 
the priority (P) of the batch job, the 
execution time (T), and central memory (CM) 
limits. 


Guaran- 
teed Cost/ 
Response SBU* 
Priority Meaning Time, hr dollars 
P1 Weekend <48 0.01 
run rr 
P2 Overnight <24+ 0.01 
run 
P3 Slower than 
default <H+ 0.04 
PY Default <2+ 0.12 
P6 Faster than 
default <1/2+ 0.13 


* Current Control Data Corporation (CDC) 
costs (FY 85). 
+ Excluding execution time. 


T: SBU account block limit. 
This is the job execution time limit (in 
seconds). 


CM: Central memory field length limit is 
the maximum amount of memory you can use 
for your batch job (in octal words). 


This form of the user card specifies that 
the user identification number being used 
in the current session will be assigned to 
this batch job. 


The charge card identifies the charge code 
being used for the current session and 
bills that charge account for the cost of 
the batch job. 


The following "GET" commands instruct the 
computer to access certain permanent 
files. Using this form of the "GET" 
command assumes that all files, except 
RCPWAVE, reside in the file space assigned 
to the user number being used in the 
current session. 


D2 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 
(13) 
(14) 


(CUS) 


(16) 


(17) 


GET , SR=RCPWAVE/UN=CER@Q2. 


GET , UPDFL=RCPUPDT. 


GET, TAPE7=RCPDATA. 


GET , TAPE8=RCPDEPT. 


UPDATE(I=SR,N-NEWLIB,F,L-0) 


UPDATE(P=NEWLIB, I=UPDFL, 
C=NEWPLUS ,L=0),F) . 


FTN5, I=NEWPLUS, B=BIN,L=OUTPUT. 


LOAD, BIN. 
EXECUTE. 
REWIND, TAPE6 , OUTPUT. 


REPLACE, TAPE6=RCPPRNT. 


REPLACE , OUTPUT=RCPOTPT. 


COST. 


The following statements instruct the 


computer to: 

Get the permanent file RCPWAVE from 

user number CER@Q2 file space and give 
it the local file name SR for this batch 
job. 

Get the permanent file RCPUPDT and 

give it the local file name UPDFL for 
this batch job. 


Get the permanent file RCPDATA and 
give iit the local file name: TAPE/ for 
this batch job. 


Get the permanent file RCPDEPT and 
give it the local file name TAPE8 for 
this batch job. 


Take the local file SR and create a full 
(F), new program library NEWLIB from it, 
but do not list (L=@) the file NEWLIB. 


Take the program library NEWLIB and add 
to it the update information contained 
in file sUPDELe “Doynot- list, the ruil 
combined program library called NEWPLUS. 


Compile the file NEWPLUS. Save the 
relocatable binary code in a local file 
BIN and the compiled listing in a local 
file OUTPUT. 


Load the local file BIN. 
Run the program BIN. 


Rewind the local files TAPE6 and OUTPUT 
which generated during compilation and 
execution of the program. 


Take the information from local file 
TAPE6 and save it in the permanent file 
RCPPRNT. If a permanent file RCPPRNT 
already exists, it will be replaced with 
this new information. 


Take the information from local file 
OUTPUT and save it in the permanent file 
RCPOTPT. If a permanent file RCPOTPT 
already exists, it will be replaced with 
this new information. 


The following statements instruct the 
computer to: 

Estimate the cost of this batch job. 
This cost excludes the CDC discount to 
US Army Corps of Engineers users. 


D3 


(18) 


(19) 


(20) 


DAYFILE,L=RCPDAYF. 


REPLACE, RCPDAYF. 


EXaary 


Create a diary of events for this batch 
job and give the diary a local file name 
RCPDAYE. 


Save RCPDAYF as a permanent file. Ifa 
permanent file RCPDAYF already exists, 
it will be replaced with this new 
information. 


Stop processing information if no errors 
have occurred. If errors have occurred 
prior to reaching line 20, the system 
jumps to this line then continues 
downward in what is called "error 
processing." Processing beyond this 
point consists of saving output and 
diary files as was done in lines 14 
through 19. 
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APPENDIX E: LINE-BY-LINE DESCRIPTION OF JOB CONTROL LANGUAGE 
FOR EXECUTING RCPWAVE ON THE CYBER 205 COMPUTER 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


/ JOB 


JOB, P3,T1000,CM377000. 


/USER 


/CHARGE 


**Tmportant** 


GET , SR=RCPWAVE/UN=CER@Q2. 


GET ,UPDFL=RCPPUPDT. 


UPDATE, 1=SR,N=NEWLIB,F,L=0. 


UPDATE, P=NEWLIB, I=UPDEL, 
C=NEWPLUS,L=0,F. 


REPLACE ,NEWPLUS. 


COPYBR, INPUT, SENDJOB, 1. 


The job card tells the computer that this 
is a batch job. 


The resource card informs the computer of 
the priority (P) of your batch job, and 
the execution time (T) and central memory 
(CM) limits. 


This form of the user card specifies that 

the user identification number being used 

in the current session will be assigned to 
this batch job. 


The charge card identifies the charge code 
being used for the current session and 
bills that charge account for the cost of 
the batch job. 


The following "GET" commands instruct the 
computer to access certain permanent 
files. Using this form of the "GET" 
command assumes that all files, except 
RCPWAVE, reside in the file space assigned 
to the user number being used in the 
current session. 


The following statements instruct the 
computer to: 


Get the permanent file RCPWAVE from user 
number CER@Q2 file space and give it the 
local file name SR for this batch job. 


Get the permanent file RCPUPDT and give it 
the local file name UPDFL for this batch 
job. 


Take the local file SR and create a full 
(F) new program library NEWLIB from it. 
Do not list (L=0) the file NEWLIB. 


Take the program library NEWLIB and add in 
the update information contained in file 
UPDFL. Do not list the full combined 
program library NEWPLUS. 


Save NEWPLUS as a permanent file. Ifa 
permanent file NEWPLUS already exists, it 
Will be replaced with this new 
information. 


Copy one record from the local file INPUT* 
to the local file SENDJOB. 


*INPUT begins after the EOR card 

(card 19). One record is defined as all 
the cards in INPUT until a /EOR or /EOF is 
encountered (Card 35). 


E2 


(13) 


(14) 


(15) 


(16) 


(19) 


(20) 


(250)) 


(22) 


SUBMIT, SENDJOB,T. 


DAYFILE,L=RCPDAYF. 


REPLACE, RCPDAYF. 


EXIT. 


/EOR 
JOB205. 


USER(U=<205 USERNAME>, 
PA=<205 PASSWORD>)ADY. 


RESOURCE( JCAT=P3,TL=2000). 


The following statements instruct the 
computer to: 


Give the CYBER 205 the batch of commands 
in file SENDJOB. 


Create a diary of events for this batch 
job and give the diary a local file name 
RCPDAYF. 


Save RCPDAYF as a permanent file. Ifa 
permanent file RCPDAYF already exists, it 
will be replaced with this new 
information. 


Stop processing information if no errors 
have occurred. If errors occur prior to 
reaching the EXIT card, the system jumps 
to this line then continues downward in 
what is called "error processing." 
Processing beyond this point consists of 
saving the diary file as was done in lines 
14 and 15. 


End of record card. All cards following 
ecard 19 are part of local file INPUT. 


The job card tells the computer that this 
is a batch job. 


The user card identifies the CYBER 205 
user for billing purposes. 


The resource card informs the computer of 
the priority (P) and time (TL) of the 
batch job. 


Guaran- 
teed Cost 
Response SBU* 
Priority Meaning Time, hr dollars 
P2 Overnight <244+ 0.08 
run 
P3 Slower than <H+ OF 2 
default 
Py Default <2+ 0.15 
P6 Faster than <1/2+ 0.18 
default i 


* Current Control Data Corporation (CDC) 
costs (FY 85). 
+ Excluding execution time. 
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(23) 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 
(30) 


(31) 


(33) 


/CHARGE 


LINK, GET ,NEWPLUS/DD=C6, 
UN=<CYBER 865 USERNAME> 
,FM=KOE, AF=AF205. 


LINK, GET, TAPE7=RCPDATA/DD=C6 


UN=<CYBER 865 USERNAME> 
,FMd=KOE, AF=AF205. 


LINK, GET , TAPE8=RCPDEPT/DD= 
C6,UN=<CYBER 865 USERNAME> 
,FM=KOE, AF=AF205. 


FORTRAN, I=NEWPLUS. 


LOAD ,ON=GO/6000,L=M205, 
GRLPALL=. 


GO. 


LINK , REPLACE(TAPE6=RCPPRNT/ 
UN=<CYBER 865 USERNAME > 
, FM=KOE, DD=C6, AF=AF205) 


SUMMARY. 


VEOK 


T: SBU account block limit. This is the 
job execution time limit (in seconds). 


The charge card identifies the charge 
being used for the current session and 
bills that charge account for the batch 
job. 


The following statements instruct the 
computer to: 


Get the permanent file NEWPLUS from the 
CYBER 865 user account. AF205 is your 
direct access file containing validation 
information. 


Get the permanent file RCPDATA from the 
CYBER 865 user account and give it the 
local file name TAPE7. 


Get the permanent file RCPDEPT from the 
CYBER 865 user account and give it the 
local file name TAPES. 


Compile the file NEWPLUS. The relocatable 
LGO by default, and the compiled listing 
is on local file OUTPUT by default. 


Load the relocatable binary code. 


Run the program. 


Take the information on TAPE6 and save it 
in the permanent file RCPPRNT. Ifa 
permanent file RCPPRNT already exists, it 
will be replaced with this new 
information. 


The summary card tells the computer to 
list out all system usage (i.e., blocks, 
disks, SBU's) for this batch job. 


End of file marker. 
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APPENDIX F: RCPWAVE PROGRAM LISTING 


“1TO ON ees Ps 


PROGRAM RCPWAVE( INPUT, OUTPUT, TAPE7, TAPEG, TAPES, TAFE3) HAIN 
CRAKAAKKARAARARAKARARKAKARKARKARAARARARRARRRARRARRARRARARRARRARARRARKAAK MAIN 


C MAIN 
Ckk NOTE kk ‘TQ’ AND ‘JQ’ DEFINE THE LENGTH OF THE ONE AND TWO HAIN 
C DIMENSIONAL ARRAYS USED IN THE PROGRAM. ‘IQ’ CORRESPONDS HAIN 
C TO THE X DIRECTION (ON-OFFSHORE) ANI ‘JQ’ 19 THE Y MAIN 
C DIRECTION (LONGSHORE : HAIN 
C MAIN 
CARKARARAAAAARKARRAARARARARARARRARARARARRARARARARARARARARARRARARARRRARAR MAIN 
C MAIN 
C PARAM 
PARAMETER ( 1Q=95, JQ=95) PARAM 
C PARAA 
C MAIN 
Ckk NOTE *k THE PRODUCT [0 & J@ MUST BE LESS THAN 9025 IF THE AODEL HAIN 
C AFPLICATION IS TO BE RUN ON THE FRONT END CYBER 863 HAIN 
C MACHINE. ANYTHING LARGER AUST BE RUN ON THE CYBER 205. MAIN 
C HAIN 
CRARAAAKARKARAKKARKARAAKARARRARRARRARARRARRARRAARARARRARKAARRARKARRARKAKK MAIN 
C HAIN 
Ckk DEFINITION OF IMPORTANT VARIABLE ARRAYS USED THROUGHOUT THE PROGRAM HAIN 
C MAIN 
C 2 - WAVE ANGLE HAIN 
C SiC =U SINGZ) HAIN 
C CO - C05(Z) HAIN 
C H - FUNCTION OF THE WAVE AMPLITULE MAIN 
C CCG - PRODUCT OF THE WAVE CELERITY AND THE GROUP VELOCITY = HAIN 
C D - TOTAL WATER DEPTH RELATIVE TO SOME DATUM MAIN 
C RKA - WAYE NUMBER DEFINED BY THE DISPERSION RELATION MAIN 
C GRDK - GRADIENT OF THE WAVE PHASE FUNCTION MAIN 
C XMUC - SCALE FACTOR RELATING REAL SPACE X GRID DISTANCES MAIN 
C TO MAPPED SPACE X GRID DISTANCES AND DEFINED AT THE MAIN 
C GRID CENTER MAIN 
C XMUS - SCALE FACTOR RELATING REAL SPACE X GRID DISTANCES MAIN 
C TO MAPPED SPACE X GRID DISTANCES ANT DEFINED AT THE AAIN 
C GRID SIDES MAIN 
C YMUC - SCALE FACTOR RELATING REAL SPACE y GRID [DISTANCES MAIN 
C TO MAPPED SPACE i GRID DISTANCES ANI DEFINED AT THE HAIN 
C GRID CENTER MAIN 
C YMUS - SCALE FACTOR RELATING REAL SPACE Y GRID [ISTANCES MAIN 
c TO MAPPED SPACE Y GRID DISTANCES AND DEFINED AT THE AAIN 
C GRID SIDES MAIN 
C XX - REAL SPACE X GRID DISTANCES MEASURE! FROM THE GRID SAIN 
C ORIGIN 10 THE GRID CENTER MAIN 
C {Y - REAL SPACE Y GRID DISTANCES MEASURED FROM THE GRID AAIN 
C ORIGIN TO THE GRID CENTER MAIN 
( HAIN 
CAAAKARAKKAARAKARAKAARAARARARKARARRRRERARRARARARRARARARARRARARAARRARARKK HAIN 
C HAIN 
Ck DEFINITIONS OF INPUT DATA TO BE READ INTO THE PROGRAM xx MAIN 
C MAIN 
C ITITLE - ARRAY CONTAINING THE APFLICATION TITLE (70 MAIN 
c CHARACTERS OR LESS) MAIN 
C H - NUMBER OF GRID CELLS IN THE X DIRECTION HAIN 
C DX - GRID SIZE IN THE X DIRECTION MAIN 
C N - NUMBER OF GRID CELLS IN THE Y DIRECTION MAIN 
C DY - GRID SIZE IN THE Y DIRECTION MAIN 
C kKNOTEAX G - GRAVITATIONAL CONSTANT WHICH DETERMINES THE UNITS OF MAIN 


372) 


Ln 4 oo bo 


(SS see oe ao 
2 QO On ee Oo me Hrd ODO |r Ss 


— 
ao 


tae 
an 


to 
= 


t2 t2 
co bro 


24 


Coe Jitee). (20 


spktsa el clalsr) (re) Creditor (sel Cel tel ese let se) Csr) (se) 


(se) (oe) tr) ose)-(se J toy) Goedel) is) 


Gz} 


C 


CNTRANG - 


NCASES 


DLEVEL 


) 


HDEEP - 
IDEEP - 
ZDEEP - 
LY 
IP2 
INC 
JFL = 
IPB 


‘ 


IGKID - 


AKNOTEAR 


IREF 


COMMON/ANGLES/ 
COMMON/WVH/H(T 
COMMON/DEPTHS/ 
COMHON/WAVNUM/ 
COMHON/ CONST,’ 
AMN,M1,MZ NAL, 
COMMON, FRINTC/ 
COMNON/CONS3/T 


ALL THE INPUT AND OUTPUT VARIABLES 

APPROXIMATE ANGLE THE OFFSHORE CONTOURS MAKE WITH 

THE Y AXIS (AIDS IN A BETTER INITIAL GUESS OF THE 
SOLUTION COMPUTED USING SNELL’S LAW) 

NUMBER OF INDIVIDUAL DEEP WATER HEIGHT-PERIOD- 
DIRECTION CONDITIONS CONSIDERED 

CONSTANT WATER LEVEL ADDED TO OK SUBTRACTED FROM 

THE ENTIRE DEFTH MATRIX 

DEEP WATER WAVE HEIGHT INPUT CONDITIONS 

DEEP WATER WAVE PERIOD INPUT CONDITIONS 

DEEP WATER WAVE DIRECTION INPUT CONDITIONS 

STARTING VALUE OF I FOR THE PRINTED OUTPUT 

ENDING YALUE OF I FOR THE PRINTED OUTPUT 

INCREMENT OF I FOR THE PRINTED OUTPUT 

STARTING VALUE OF J FOR THE PRINTED QUTPUT 

ENDING VALUE OF J FOR THE PRINTED OUTPUT 

(THE J INCREMENT 15 ALWAYS 1) 

=0 IF CONSIDERING A CONSTANT SIZED RECTANGULAR GRID 
=] IF CONSIDERING A YARIABLY SIZED RECTANGULAR ‘3k 1D 
GENERATED USING THE MAPIT PROCEDURE IN CMSGRID. 
FOR ADDITIONAL INFORMATION ABOUT THE YARIABLY SIZED 
GRID INPUT SEE THE COMMENTS IN ‘SUBROUTINE GRID’ 
=] DIFFRACTIVE EFFECTS INCLUDED 

=Q DIFFRACTIVE EFFCTS IGNORED(PURE REFRACTION) 


KAAKKARKKARAARARAKARARARAARARARRARAAKARRARRERARAARERRRARARAARARRERARRRR 


Z(10,JQ),S1( 18, JQ) COC 10, JQ) 

Q,J@) ,CCG( 10. JQ) 

DC 1G, JQ) 

RKA(1Q, JQ) ,GRDK (10. JQ) 
G,P1,P12,RAL,HCONVR,SCONVR DX, DY ,DX2,DY2,7,OMEG, 
NM2, IWVET, IDRY, IWETP1 CNTRANG,HFACT DLEVEL 

IP], 1P2, INC, JP1,JP2,DMULT,HMULT, ZMULT, RKMULT 
GRID, IREF, [TAHX, ITHMX. IDIFF 


COMMON/MUSC/XMUC (IQ) , XMUS( 10), YMUC( JQ) , YMUS( JQ) 


COMMON/COOR/XX 
COMHON /TRANE/ 
DIMENSION HDEE 


(TQ) YY (JQ) 
DECAY ,STABL, IBRK( JQ) , IRRKM( JQ) 
P(200) , TDEEP(200) ,2DEEP (200) .LTITLE(9) 


Ckkkkkakh READ INPUT DATA FROM FILE CODE ” 


C 


C 


READ(? 601) (LT 
601 FORMAT (9A8) 


ITLE(L),L=1,9) 


READ (7, 100)M,DX,N,DY,G,CNTRANG,NCASES , DLEVEL 


100 FORMAT(1S,F10. 


4,15,F10.4,2F10.2,15,F10.2) 


DO 602 L=1,NCASES 
READ(7,103)HDEEP(L) , TDEEP(L) ,ZDEEP(L) 


103 FORMAT(3F10.2) 

602 CONTINUE 
READ(7, 101) IFL 

101 FORMAT(SI5) 


»IPd, INC, JP1,JF2 


READ(?, 102) IREF, IGRID 


102 FORMAT(215) 


CkkkkAAKK WRITE OUT INPUT DATA ON FILE CODE 6 


C 


WRITE(6,15) 


15 FORMAT(/////,2X,/REG TONAL COASTAL’,3X, 
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113 
119 
120 
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PROCESSES WAVE TRANSEORMATION, 
43X,/M ODE L,/,43X,/0 ROPWAVE)’,////) 
WRITE(6,609) (LIITLE(L) ,L=1.9) 
609 EORMAT(25X,9A8, ///) 
WRITE(6,55) 
55 FORMAT(,’,1X,/MOUEL INPUT:’,//) 
IE‘ IGRID.GT.0)G0 TO 300 
WRITE(6, 41) 
4] FORMAT(1X, ‘UNIFORMLY SIZED, RECTILINEAR GRID MESH’.,’) 
80 10 301 
300 WRITE(6, 42) 
42 FORMAT(1X,VARIABLY SIZED, RECTILINEAR GRID MESH’ , /) 
301 WRITE(6,5)4,0X,N,DY 
5 FORMAT(5X, ‘Y - DIRECTION (ON-OFESHORE)’,5X,1S,’ CELLS, EACH’, 
KE10.2,’ IN LENGTH’./,5X,’Y - DIRECTION (LONGSHOKE) ’,5X, 15. 
k’ CELLS, EACH’,E10.2,° IN LENGTH’, /) 
WRITE‘, 43)6 
43 FORMAT(1X, THE ACCELERATION OF GRAVITY IS’,f8.2,’. THESE UNITS 
ADETERMINE THE UNITS OF ALL INPUT AND OUTPUT VARIABLES ’,,) 
DO 604 L=1,NCASES 
WRITE(G,44)L yHUEEP(L) , 1DEEP(L) , 2DEEP(L) 
604 CONTINUE 
44 EORMAT(’ THE DEEP WATER WAVE PARAMETERS FOR CASE’.13,’ ARE:’, 
45X, ‘HEIGHT=’ ,£7.3,5X, ’PERIOD=‘,£7.3,5%, /ANGLE=’ ,£8.3) 
WRITE (6,45) CNTRANG 
45 FORMAT(/,1X, ‘THE OFFSHORE CONTOURS MAKE AN ANGLE OF’ ,E7.2, 
k’ DEGREES WITH THE Y AxI3’./) 
WRITE (6, 193 


19 FORMAT(1X, THE BATHYMETRY MATRIX IS VARIABLE IN BOTH HORIZONTAL’ 


k,‘ DIRECTIONS ANI! WAS READ FROM FILE CODE 3’,/) 
WRITE(5,20)DLEVEL 
20 FORMAT(1X,‘A WATER LEVEL CHANGE OF’,F7.2,’ WAS ADDED TO THE’. 
A’ ENTIRE BATHYMETRY MATRIX’) 
C 
Chrkhrkhh CONSTANTS USED IN THE PROGRAM 
C 
P1=3.1415927 
FIZ=PI &2.0 
RaAD=180.0/P1 
CNTRANG=CNTRANG/ RAD 
Ha=M-2 
H1=H-1 
NM1=N-1 
NHQ=N-2 
DXE=DXKL. 
DY2=DYAc. 
WIA=1.0 
WIH=1.0 
TAU=0.167 
BETA=0.167 
IDRY=1 
IWET=2 
IWETP1=IWET+1 
C 
Ckkkkkkkk ‘STABL’ AND ‘DECAY’ ARE USED IN THE WAVE BREAKING 
CrkrkkakK ROUTINE 
C 
STABL=0.4 
DECAY=0.2 
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C 
ChkRARAAR HCONVR’ AND ‘SCONVR’ ARE THE CONVERGENCE CRITERIA USED IN 
ChAkARAAA THE HEIGHT AND ANGLE ITERATIVE SOLUTION SCHEMES 
L 
HCONVR=0..0005 
SCONVR=0.90025 
IF(G.G7.9.7.AND.G.LT.9.9)HCONVR=HCONVRAO. 3048 
1F(G.G7.979.0, AND. G.LT.990.0)HCONVR=HCONVRA3O. 48 
C 


Cakkakkah ‘ITHMX’ AND ‘TTAMX’ CONTROL THE MAXIMUM NUMBER UF ITERATIONS 


CkakARKKK ALLOWED IN THE WAVE HEIGHT AND ANGLE SOLUTION SCHEMES 
Chkkhakak “IDIEF’ CONTROLS THE NUMBER OF ITERATIONS ALLOWED IN THE 
Craarctah ITERATIVE DIFFRACTION SCHEME 
C 

TTHMX=5¢ 

TTAMX=50 

INIFF=15 


st 


L 


CakakkAAA THE FOLLOWING MULTIPLICATION FACTORS CONTROL THE ACCURACY 
CAkkAARtA OF THE PRINTED QUIPUT‘NOT THE ACTUAL COMPUTATION). 


CRAARARRK DEPTH (1) ~ DMULT 
CARRAKARR WAVE HEIGHT (H) - HMULT 
CARARAARR WAVE ANGLE (2) - ZMULT 


ChARRAKKK WAVE NUMBER (GRADK) - RKMULT 
CkkkARKAK THE VARIABLES ARE FIRST MULTIPLIED BY THE SCALE FACTORS 
Crkthkakk KELOW, THEN PRINTED QUT IN INTEGER FORM 
L 
UMULT=1.9 
HAULT=10.6 
DMULT=1.9 
REMULT=1000.0 
C 
CAkkkaAAR DEFINE THE GRID SCALE FACTORS AND DISTANCES FOR A CONSTANT 
Charakhhe SIZED RECTANGULAR GRID 
C 
10 2 I=1,4 
XMUC(1)=1.0 
XAUS(1)=1.0 
XX(T)=DX4(1-0.5) 
2 CONTINUE 
0 3 J=1.¥ 
(AUC(I)=1.9 
YMUStJ)=1.0 
WD) =0fACI-9.5) 
3 CONTINUE 


4 


ry 


i 


Sy 


Chkkkhkak CALL SUBROUITNE GRIL TO READ IN THE VARIABLE GRID MAPPING 
Chkkkkkkk INFORMATION FROM FILE CODE 3 (X DIRECTION FIRST, THEN 1) AND 
Chakthhkh GENERATE THE VARIARLE GRID SCALE FACTORS AND [DISTANCES 
C 

TE(IGRID.EQ.1)CALL GRID 


U 


Ckkkkakkh CALL SUBROUTINE DEPTH TO GENERATE THE WATER DEPTH MATKIX 
C 

CALL DEPTH 
C 
Chkkkkkhkh ADD [LEVEL TO ALL THE WATER DEPTHS AND THEN REQUIKE 
Chkkkktkk ALL DEPTHS BE GREATER THAN THE MINIMUM DEPTH ‘DMIN’ 
Ckkkkkake WHERE “DMIN’ IS 1.0 FEET OR A METRIC EQUIVALENT 
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DMIN=1.6 
1E(G.G7.9.7.AND..LT.9.9) DMIN=0. 2048 
1E(G.67.970.0.AND.G.LT.990.0)DHIN=30. 48 
DO 87 I=1,4 
10 87 J=1,N 
D(I,J)=0(1, J} #DLEVEL 
TE(D(1,J).LT.DMIN)D( 1, J)=DMIN 

87 CONTINUE 

G. 

Chkkkkkth SET THE BOUNDARY CONDITIONS FOR THE WATER DEPTHS 

C 
DO 86 I=1,4 
D(I,1)=0(1,2) 

D4 1.N)=D(1,NAL) 
86 CONTINUE 
CALL ONBC(D) 


c 
CkkakkkaR PRINT THE TOTAL WATER DEPTH MATRIX ‘ALL DEPTHS MULTIPLIED BY 0 
C 

CALL POUT(IP],IP2,INC,JP1.JP2,’ WATER DEPTHS ",D,DAULT) 
C 


Ckkkkkkhk ITERATE OVER THE NUMBER OF WAVE CONDITIONS CONSIDEREL 
C 
DO 605 L=1,NCASES 
HO=HDEEP(L) 
T=TDEEP(L) 
OMEG=P 12/T 
HEACT=2. 0XOHEG/G 
C 
Cakkkkakk ‘HEACT’ I5 A FACTOR T9 CONVERT FROM WAVE HEIGHT TO THE 
CrkkAAKAA AMPLITUDE FUNCTION ‘H’ 
C 
A=ZDEEP(L) 
WRITE(6,606)L 
606 FORMAT(////",° WAVE CONDITION’,1,/) 
WRITE(6,44)L,HO,T,A 
A=A+180.0 
0 507 J=1,N 
IBRK(J)=0 
IBRKH(J)=0 
607 CONTINUE 
CALL REFDIE(A,HO,WIA.TAU.WIH, RETA) 
605 CONTINUE 


99 STOF 
END 


CARAAAKAARAKAKRARARRAARARARAARARKARARRARRAARRRAARRARAR AR ARRARRRARRARRA RA 


C 
SUBROUTINE REFDIE(THETAO,HH,WIA, TAU,WIH, BETA? 
C 


CRRARARKARAAKEARARKARAAARARRARARARARARERARARRARRARRARRARARAARRRRRARRRR RA 


C 

C THIS SUBROUTINE CONTROLS THE ITERATIVE SOLUTION SCHEME AND 
C PRINTING OF THE WAVE INFORMATION 

C 


CARKAAKAARARARARARARAARAAARARKARRARRARRARRARRARRARRARARARRRRAR AR ARARARRR 


C 
PARAHETER( I0=95,, JQ=95) 
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COMMON/ANGLES/2( 10. JQ) ,51(10,JQ) CO’ 12, JQ) 
COMMON/WVH/H(1,JQ) CCG 10,1) 

COMMON/ DEPTHS’ 10, JQ: 

COMMON, WAVNUM/RKA: IQ, JG) GRDK( 10, JQ: 


COMMON/CONST,’ G, PI, PID, KAL,HCONYR, SCONE, DX, 0Y.0X2 02,1. OMEG, 


AM N, HL, MO,NM] NMC, IWET, IDRY, IWETFi.CNTRANG HEACT [LEVEL 
COMMON/PRINTC,/IP1,IP2, INC, JP1,JP2,DMULT HAULT.ZAULT, KKMULT 
COMHON/CONS2/ IGRID, IREF, ITAMY, ITHAX, IDIFE 
COMMON /TRANE, DECAY, STABL, [BRK(JQ), IBRKK( JQ) 

DIMENSION DUMi‘ 10, JQ) ,GREOLD JQ: 
C 
CaakARRAA CALL SUBROUTINE SNELL TO OBTAIN AN INITIAL GUESS GF THE 
CkkaARARK WAVE HEIGHTS AND ANGLES USING SNELLS LAW 
C 

CALL SNELL |‘ THETAG,’RAD) HH) 

CALL ONKC(H) 

CALL ONBC (2) 

CALL ONBC(CO: 

CALL ONBC{51} 

CALL ONBC(RKA) 

CALL ONBC(SRDK! 

CALL ONBC{ CCI) 
C 


Crkakkakh REMOVE THE “0 TO i232 STATEMENT TO PRINT OUT THE SNELL’S LAW 


Carrkkekr SOLUTION 
C 
40 TO 125 
177 CONTINUE 
00 13 I=1,4 
110 13 J=1.N 
QUHG=2(1,J)-PI 
DUM (1,J)=DUN3*RAD 
12 CONTINUE 
C 
Ckakaktak ANGLES HEIGHTS, AND WAVE NUMBERS AXE AULIIPLIED Sy ZMULT. 


CrkARKAAA HMULT, AND REMULT RESPECTIVELY IN THE FRINT OUT. THESE MUST BE 
Ckkkkkakk CHANGED INTERNALLY WITHIN THE MAIN FROGRAM Ir OTHER ‘ALUES 


Cakakrare ARE DESIRED 


CALL POUT(IFL,IPS,INC,JPL,JPo,’ WAVE ANGLES = (DEG »DUHL, 


«ZMULT: 

1G 79 I=1,é 

10 73 J=1.N 

DUAL 1.5)=H¢ 1,3) xHEACT 
79 CONTINUE 


CALL POUT(IP1,IP2,INC,JP1,JP2,’ WAVE HEIGHTS * DUM1, 
1HAULT} 
CALL POUT(IP1,IP2,INC,JF1,JFo,’ WAVE NUMBER ' RRA, 
ARKMULT) 
122 CONTINUE 
C 
CrkkkaARK ROW BY ROW MARCHING LOOP 
c 
DO 10 IL=H2, IWETP1,-1 
ILH1=IL-1 
C 
Ckkkkkakk COMPUTE THE WAVE ANGLES ALONG A ROW 
C 
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CALL ANGLE(IL, IL, ITAHX,WIA, TAU, SCONVR) 


C 
Chkkkkakkk COMPUTE THE WAVE HEIGHTS ALONG A ROW 
C 

CALL HEIGHT (IL, IL, ITHAX,WIH, BETA,HCONVR) 
C 


TE( IREF.EQ.0) GO 10 114 
DO 119 J=1,N 
GRDOLD(J)=GRDK( ILH1, J) 
119 CONTINUE 
DO 110 LLL=1, IDIEF 
C 
Ckrkkkkkk COMPUTE NEW GRADIENT OF THE WAVE PHASE FUNCTION ALONG A ROW 
C 
CALL GRADK( ILH1, ILM1) 
EPSK=0.95 
DO 112 J=1,N 
GRDK ( ILM], J)=EPSKAGRDK ( ILM), J)+(1.0-EPSK) AGRDOLD( J) 
112 CONTINUE 
TFLAG=1 
HO 111 J=1,N 
TE(ABS(GRDK ( ILM1,J)-GRIGLD(J)) .GT.0.00254ABS(GRDK ( ILK1,J))) 
kIFLAG=0 
111 CONTINUE 
{0 113 J=1,N 
GRDOLD (J) =GRDK( ILM1,J) 
113 CONTINUE 
C 
Cakkkakkk RECOMPUTE THE WAVE ANGLES ALONG A ROW 
C 
CALL ANGLE(IL, IL, ITAMX,WIA, (TAUAL.0) ,SCONVR) 


C 
Ckkkakkth RECOMPUTE THE WAVE HEIGHTS ALONG A ROW 
C 
CALL HEIGHT( IL, IL, ITHMX,WIH, (BETAX1.0) ,HCONVR) 
TE(TELAG.EQ.1) GO TO 114 
110 CONTINUE 
WRITE(6,117) ILA] 
117 FORMAT(1X, ‘CONVERGENCE TOWARD A SOLUTION FAILED ON ROW’, I3) 
114 CONTINUE 


C 

Crkkrrkkh UPDATE BREAKER INDEX 

C 

BO 115 J=1,N 
IBRKM(J)=IBRK( J) 
DUM1( ILH1,J)=IBRK(J) 
IBRK(J)=0 
115 CONTINUE 
C 


10 CONTINUE 
CALL ONBC(GRDK) 
CALL ONBC(H) 
CALL ONBC(Z) 
CALL ONBC(CO) 
CALL ONBC(SI) 
C 
Crrkkkakk PRINT WAVE INFORMATION 
C 
Chkkkkkkk BREAKER INDEX ( O=NON-BREAKING , 11=BREAKING ) 
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DO 219 J=1,N 

DUM) (1, J)=DUML(2, J) 

DUMI(M,J)=0.0 

DUAL (M1,J)=0.0 

DUML(H2,J)=0.0 

DO 219 I=1,4 

TE(DUMI (I,J) .EQ.1.0)DUMI(I,J)=11.0 
219 CONTINUE 


CALL POUT(IP1,IP2,INC,JP1,JP2,’ BREAKER INDEX “,DUM1,1.0) 


CrxaKKARA ANGLES HEIGHTS, AND WAVE NUMBERS ARE MULTIPLIED BY ZMULT, 


CarkaARRAK HMULT, AND RKMULT RESPECTIVELY IN THE PRINT OUT. THESE MUST 
Crkkkkkkk BE CHANGED INTERNALLY WITHIN THE MAIN PROGRAM IF OTHER VALUES 


Chkkkhkhk ARE DESIRED 
c 
DO 23 I=1,H 
DO 23 J=1,N 
DUM=2(1,J)-PI 
DUM] ( 1, J)=DUMSARAD 
23 CONTINUE 
CALL POUT(IP1,IP2, INC,JP1,JP2,’ WAVE ANGLES © (DEG)’,DUML, 
AZMULT) 
C 
Ckkkkkkkh CONVERT EROM AMPLITUDE FUNCTION TO WAVE HEIGHT 
C 
DO 81 I=1,4 
DO 81 J=1,N 
DUML(I,Ji=H( I,J) AHEACT 
81 CONTINUE 
CALL POUT(IP1,IP2,INC,JP1,JP2,’ WAVE HEIGHTS ’ DUM, 
1HMULT) 
CALL POUT(IP1,IP2, INC, JP1.JP2,’  4AVE NUMBER ’ GRDK, 
ARKMULT) 
99. RETURN 
END 


CRARAAKAKARKAAKARAAAAAARARARARAAAAARAARARRAARARRARRARARAARRARRARRRAARRRR 


C 
SUBROUTINE GRADK( ISTART, TEND) 
C 


CARAKAARAKARARAARAAKARARAAARARAARRARAARARKARAARARARAARRARERARRRRAAARARRKK 


C 

Cc THIS SUBROUTINE COMPUTES THE UPDATED VALUES OF THE GRADIENT 
C OF THE WAVE PHASE FUNCTION ALONG A GIVEN ROW 

c 


CAARARARAAARAAAARARRARKARKARRARAAAAARARARAARARARRARARRARRARRRARARRARK ARK 


C 
PARAMETER ( 10=95, JQ=95) 

C 
COMMON/ANGLES/Z(1G,J@).S1(1G,JQ) ,CO( I,J) 
COMMON/DEPTHS/D( 1G, JQ) 
COHMON/WVH/H( 18, JQ) ,CCG( 18, JQ) 
COMMON/WAVNUHM/RKA( TQ, JQ) ,GRDK( 10, JQ) 
COMMON/CONST/ G,PI,PI2,RAD,HCONVR,SCONVR,DX,DY,DX2,D¥2,T,0MEG, 
AM N,H1,M2,NM1,NM2, IWET, IDRY, IWETP],CNTRANG , HEACT, DLEVEL 
COMMON/CONS2/IGRID, IREF, ITAHX, ITHMX, IDIFE 
COMMON/MUSC/XMUC( 1), XMUS( TQ) , YMUC( JQ) , YMUS (JQ) 
COMMON /TRANE/ DECAY,STABL, IBRK( JQ), IBRKM(JQ) 
DIMENSION DUM1(10,JQ) ,x(JQ) 


Eg 


REFDIE 
REFDIE 
REFDIE 
REFDIE 
REFDIF 
REFDIE 
REEDIE 
REFDIF 
KEEDIE 
REFDIE 
REFDIE 
REED IEF 
REEDIE 
REFDIEF 
REFDIF 
REFDIF 
REEDIF 
REFLIIF 
REEDIF 
REFDIE 
REEDIF 
REFDIE 
REEDIE 
REEDIF 
REED IF 
REFDIF 
REEDIF 
REEDIF 
REFDIF 
REEDIF 
REFDIE 
KEFDIEF 
REFDIE 
REFDIE 
REEDIF 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
PARAM 
PARAM 
PARAM 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 
GRADK 


“1 Oo Cn & OF 2 


me 
Be worse Ow wo 


— 
& Co 


eee 
“10 cn 


se 
woo 


ty 


to brs 
[oe od 


to 


Ckxkkkkkk COMPUTE THE DIFFRACTIVE TERMS 
C 
Crkkakkkk USE 4 POINT BACKWARDS DIFFERENCE FOR THE X CURVATURE OF ’H’ 
Chkkkkkkk USE 3 POINT BACKWARDS DIFFERENCE FOR THE % GRADIENT OF ‘H,CCG’ 
Chkkkakkk USE CENTRAL DIFFERENCES FOR THE { CURVATURE OF ‘H’ 
Chkkkakkk USE CENTRAL DIFFERENCES FOR THE Y GRADIENT OF ’H,CCG’ 
C 
DO 2 I=ISTART, TEND 
DLX2=DX2AXHUC (1) 
DLXSO=(DLX240.9) AK2 
DO S J=2,NM1 
DLY2=DY24YKUC (J) 
DLYSQ=(DLY240.5) xk2 
CCGIJ=CCG(1, J) 
HIJ=H(1,J) 
HIPI=H(1, J+1) 
HJM1=H(1,J-1) 
HIPI=H(I+1,J' 
HIPZ=H(1+2,J) 
DUM6=(-3.04H1I+4. OAHIP1-HIP2) /DLX2 
DUM4=1 2. 0AHII-5. OAHIP1 +4. OAH IP2-H( I+3, 1) )/DLXSQ 
DUM4=DUM4-0.. SADUMGA (XMUC( I+] )-XMUC(I-1) )/((XMUC( 1) kk2) KDX) 
DUA2=(-3. 0aCCGII+4. OACCG( 1+], J)-CCG( +2.) )/DLX2 
DUM7=(HJP1-HIM1) /DLY2 
DUMS=(HIP1-2.0AHII+HJM1) /DLYSO 
DUMS=DUMS-0. SADUM7&( YMUC ( J+1)-YMUC(J-1))/( (YMUC(J) &&2) ALY) 
DUM3=(CCG(1,J+1)-CCG(1,J-1))/DLY2 
DUM1 (1, J)=(DUR4+DURS+ ( DUM2ADUM6+DUM3ADUM7 ) /CCG IJ) /HIJ 
= CONTINUE 
C 
Ckkkkkkkk CHECK FOR POINT TO POINT OSCILLATIONS IN WAVE PHASE GRADIENT 
C 
10 410 I=ISTART, IEND 
DUM1(1,1)=DUM1( 1,2) 
(WMI (1.N)=DUM1 (1, NAL) 
SKU=4.¢ 
SBETA=0. 25 
10 409 J=1,N 
X(J)=DUAL(I.J) 
409 CONTINUE 
0 408 J=3.NM2 
DUMZ=ABS(X(J+1)-X(J)) 
DUAS=ABS(X(J)-X(J-1)) 
DUM4=ABS (X(J+1)-X(J-1) 40.5 
TF ((DUM2+DUM3) LT. (SMUADUM4) GO TO 408 
DUMS=X(J+1)-2.0aX(J)+X(J-1) 
DUMG=X{1+2)-2.0AX(J+1)+X(J) 
DUM7=X(J)-2.04X(J-1)+X(J-2) 
TE( (DUMSADUHG) .GE.0.0.AND. (DUMSADUM7) .GE.0.0)G0 T0 408 
DUM1 (1, J)=X(J)+SBETAA(X(J+1)-2.0AX(J)+X(J-1)) 
408 CONTINUE 
410 CONTINUE 


DO 17 I=ISTART, TEND 
DO 17 J=2,NH1 


Ckkkkkkkk DO NOT INCLUDE DIFFRACTIVE EFFECTS INSIDE OR IMMEDIATELY 


CkkkkkARK ADJACENT TO THE BREAKER ZONE 
c 
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IBTEST=IBRKA( J-1)+IBRKM(J)+IBRKM(J+1)+IBRK(J-1)+IBRK(J)+ GRADK 
IBRK(J+#1) GRADK 
IE(IBTEST.GT.0) G0 0 17 GRADK 
c GRADK 
CxarkAKKK LIMIT THE CHANGE IN THE GRADIENT OF THE WAVE PHASE FUNCTION GRADK 
Carkkakkk FROM THE DISPERSION RELATION WAVE NUMBER TO 50 PERCENT GRADK 
c GRADK 
RKOLD=RKA( I,J) GRADK 
RKADD=0. 254RKOLD GRADK 
RKARG=RKOLDA42+DUM1 (1, J) GRADK 
TE (RKARG.LE.0.0) RKARG=( (OMEGAA2/G) &42) GRADK 
RKNEW=SQRT(RKARG) GRADK 
RKD TFF=RKNEW-RKOLD GRADK 
TF(ABS(RKDIFF) .LT.(0.0025ARKOLD) G0 TO 17 GRADK 
TE(ABS(RKDIFE) .LE.RKADD)GO TO 341 GRADK 
TE(RKDIFE.LT.0.0) RKNEW=RKOLD-RKADD GRADK 
TE(REDIFE.GE.0.0) RKNEW=RKOLD+RKADD GRADK 
341 GRDK(1,J)=RKNEW GRADK 
17 CONTINUE GRADE 
C GRADK 
Cxrkkkkkk SET LATERAL BOUNDARY CONDITIONS FOR THE GRADIENT OF THE WAVE SGRADK 
Cakkakkak PHASE FUNCTION GRADK 
C GRADK 
DO 38 I=ISTART, TEND GRADK 
GRDK(1,1)=GRDK(I,2) SRADK 
GRDK(1,H)=GRBK( 1,1) GRADK 
38 CONTINUE GRADK 
RETURN GRADK 
END GRADK 
Cakkkrhrakkararaahakrkraarkrrerkhaaraarka akan kariaaakkkkkkek HEIGHT 
C HEIGHT 
SUBROUTINE HEIGHT( ISTART, IEND, ITMAX,WI,ALPHA,HCNY) HEIGHT 
C HEIGHT 
Carkarkrrrraraeraak kkhnhakankaaararakaahhakhkikkakekeakkrekkk HEIGHT 
Cc HEIGHT 
c THIS SUBROUTINE ITERATES TO SOLVE FOR THE WAVE HEIGHTS HEIGHT 
C ALONG A GIVEN ROW HEIGHT 
Cc HEIGHT 
Cakkacarknakertachinkrcrararancrgnrragrkaaaakrekaaerharkerakrrkkahe HEIGHT 
Cc PARAM 
PARAMETER ( [0=95, JQ=95) PARAK 
C PARAM 
COMMON/ANGLES/Z( 10, JG) ,$1(10,J@) ,CO(1G,J&) HEIGHT - 
COMMON/WVH/H(10,J@) ,CCG( 1G, JQ) HEIGHT 
COMMON/DEPTHS/D(1@,J@) HEIGHT 
COMMON/WAVNUM/RKA( 10, JQ) ,GRDK( 19,10) HE IGHT 
COMMON/CONST/ G,P1,P12,RAD,HCONVK,SCONVK,DX,DY,DX2,DY2,7,OMEG, HEIGHT 
xM,H,M1,42,WM1,NM2, IWET, IDRY, IWETP1, CHTRANG,HEACT, DLEVEL HEIGHT 
COMMON/CONS2/IGRID, IREF , ITAMX, ITHMX, IDIFF HEIGHT 
COMMON/MUSC/XMUC( IQ), XHUS(1Q), YHUC( JQ) , YRUS(JQ) HEIGHT 
COMMON/COOR/XX(1@) , Y¥(J@) HEIGHT 
COMMON /TRANE/ DECAY,STABL, IBRK(J@) , IBRKM(J@) HEIGHT 
DIMENSION DM8(JQ) ,SLOPE(JG) HEIGHT 
Cc HE IGHT 
Ckkdckhth SOLVE DIFFERENCED FORM OF CONSERVATION OF WAVE ACTION EQUATION HEIGHT 
Cc HEIGHT 
DO 500 I=IEND, ISTART,-1 HEIGHT 
IM=I-1 HE IGHT 
DO 450 IT=1, ITHAX HEIGHT 
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DO 200 J=2,NM1 

JP=J+1 
JM=J-1 
DML=CCG( 1M, JP) AGRDK( IM, JP) ASI( IK, JP) 
DM2=CCG( IM, JM) AGRDK( IM, JM) ASI( TH, JM) 
DAS=CCG(1, JP) AGRDK (1, JP) ASI(1, JF) 
DM4=CCG(1, JK) AGRDK( 1, JH) ASI(1, JM) 
DMG=(WTA(DM1AH( TH, JP) AA2-DMZAH( TH, JH) kA) 

= /(DY2AYMUC()))+((1.0-WT)A(DMSKH( 1, JP) AX2-DM Ak 

= H(T, JM) 4x2) /(DY2KYMUC( J) )) 


DML=(COG( 1, J) AGRDK(1,J)ACO(1, J) )ACH(T, J) AH (I,J) ) 
DM2=(CCG(1,JK)AGRDK(1,JK)ACO(1, JM) )A(H(1, JM) AH(1, JK) ) 
DM3=(CCG( 1, JP) AGRDK( 1, JP)ACO(1, JP) )ACH( 1, IP) AH( 1, IP) ) 
DM4=CCG( IM, J) AGRDK( IM, J) ACO( IM, J) 
DMS=( (DMG) A(DXAXMUS( IM) ) +ALPHAADH2+ 
e (1.0-2. OXALPHA) ADH] +ALPHAADH3) /DM4 
C 
ChxkkAAKK CONVERT AMPLITUDE FUNCTION TO WAVE HEIGHT 
C 
DM7=DMSA(HEACTAHEACT ) 
C 
Chkkkkkkk CHECK FOR WAVE HEIGHTS LESS THAN ZERO 
C 
IE(DM? .LT.0.00001) DH7=0.00001 
DM7=SQRT(DM7) 
C 
Ckkkkkakk CHECK FOR WAVE BREAKING OR WAVE TRANSFORMATION 
C 
IF(DM7 .GT.STABLAD(IH,J)) THEN 
C 
Crkakkkkk CALCULATE ANGLES AND CELLS OF INFLUENCE 
c 
THETA=Z(1M,J)-PI 
IF(THETA.GE.0.) THEN 
410.54 (XX (IM) +XX( 1) )-XX( TH) 
XQ=XX(1)-XX (1M) 
Y1=YY(J+1)-YY (1) 
Y2=0 OA(YY (J) 4Y¥ (J#1) )-¥¥ (J) 
ANG1=ATAN2(Y1,X1) 
ANG2=ATAN2(Y2,X2) 


IE(THETA.GT.ANGL) THEN 
KEY=IBRK(J+1) 
IKEY=IH 
JKEY=J+1 

ELSE IF(THETA.GT.ANG2) THEN 
KEY=IBRKM(J+1) 
IKEY=1 
IKEY=J+1 

ELSE 
KEY=IBRKM( J) 

IKEY=1 
JKEY=J 
END IE 


ELSE 
X3=XX(1)-XX( 1M) 
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C 


X4=0 GR (XX (THK) +XX( 1) )-XX (1M) 
Y3=0..5A(YY (J) +¥¥(J-1))-1Y (J) 
Y4=Y¥(J-1)-YY(J) 
ANG3=ATAN2 (13, X3) 
ANG4=ATAN2(Y4,X4) 


TE(THETA.GT.ANG3) THEN 
KEY=IBRKM( J) 
IKEY=] 
IKEY=J 

ELSE IE(THETA.GT.ANG4) THEN 
KEY=IBRKM(J-1) 
IKEY=1I 
JKEY=J-1 

ELSE 
KEY=IBRK(J-1) 
IKEY=I¥ 
JKEY=J-1 

END IE 


ENDIF 


CrkkkkkKR COMPUTE INCIPIENT BREAKING WAVE HEIGHT USING METHOD FROM 
CrrkkkAAK THE SHORE PROTECTION MANUAL (WEGGEL’S EMPIRICAL METHOD) 


C 


C 


Jd=(J+(IKEY-J-1)/2) 
DIST=SQRT( (FLOAT ( IKEY- IM) ADXAXMUS( IM) ) AAZ+ 


(FLOAT (JKEY-J) ADYAYMUS (JJ) ) x42) 


SLOPE (J)=(D( IKEY , JKEY)-D( 1M,J))/DIST 
SLOPE (J)=AMAX1(0.0,SLOPE(J)) 

HS IG=1. OADH7 
Al=43.754(1.0-EXP(-19.0xSLOPE(J)) ) 
Bl=1.96/(1.0+EXP(-19. S4SLOPE(]) )) 
DEPBRK=HS1G/(B1-(A1 AHS 1G/(GATAT) ) ) 
IF(DEPBRK.GT.D(IK,J)) KEY=1 


CkARKKAKA CHECK WHETHER TO ACTUALLY TRANSFORM WAVE 


C 


C 


C 


C 


IF(KEY.EQ.1) THEN 


IBRK(J)=1 

DM1=0.3k(D( IKEY, JKEY)+D( 1H, J)) 

DM2=CCG( IKEY , JKEY )AGRDK ( TREY, JREY) AC H( IKEY , JREY ) Ax2- 
(STABLAD( IKEY , JKEY)/HEACT) Ak2)-CCG( TH, J) AGRDK( IH, J)’ 
(STABLAD( IM, J)/HEACT) Axo 

DMS=DXAXHUS (1M) ADECAY,/(2. OADM DMA ) 

DMS=(DMS+DMSADH2)/(1.0-CCG( IM, J) AGRDK( 1K, J) ADM3) 

DM7=DMSx (HEACTAHEACT ) 

TE(DM7 LT. ((STABLAD( IH, J) )Ak2)) DN7=(STABLAD( IM, J) ) 4x2 

DM7=SQRT (D7) 


ENDIE 


DAG (1) =D47 


Ckkkakirr CHECK FOR WAVE HEIGHT CONVERGENCE 


C 


DO 400 J=2,NM1 
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047=048 (J) 
HH2=H( IM, J) AHEACT 


C 

Crkkkkkkk UPDATE THE WAVE HEIGHT 

C 
EPSZ=0.9 
DA7=EPSZADM7+(1.0-EPSZ) kHH2 
SUM=SUMK+ABS ( DA7-HH2 ) 

C 


Carkaxkkk CONVERT HEIGHT BACK TO AMPLITUDE FUNCTION 
C 
H(IM,J)=DH7/HEACT 
400 CONTINUE 
C 
Crkkkkkkk SET LATERAL BOUNDARY CONDITIONS FOR WAVE HEIGHT 
C 
H(TM,1)=H(IM,2) 
H(IM,N)=H( IK, NM1) 
C 
Crkraakkk TEST FOR CONVERGENCE 
C 
TE(SUM.LT. (ELOAT(NH2) AHCNV)) GO TO 500 
450 CONTINUE 
¢ 
IF(IREF.EQ.1)G0 TO 500 
WRITE(6,475) IM 


473 FORMAT(1X, ‘RELAXATION FOR WAVE HEIGHTS FAILED ON ROW =’, 14) 


C 
500 CONTINUE 
RETURN 
END 
pvyerevercecyecevirusccrcercvryeceryersrrerververervrrerercerecceryr tse! 
c 
SUBROUTINE SNELL (THETAO,HH) 
c 
errrcoceererierrreccrrrrrererrerirrrrcrrrerirrrrerert rete rr rere rites 
C 
C THIS SUBROUTINE COMPUTES THE SNELLS LAW SOLUTION OVER THE ENTIRE 
C GRID 
C 
Ciacci kod dadicdddicdacaic oda ddd dk 
C 
PARAMETER( 10=95, J0=95) 
C 
COMMON/ANGLES/Z( 10, JQ) ,$1(10,J0) ,CO( 10, JQ) 
COMMON/WVH/H( 10, JQ) ,CCG( 10, JQ) 
COMMON/DEPTHS/D\ 10, J0) 
COMMON/WAVNUM/RKA( 10, JQ) ,GRDK( 12, JQ) 
COMMON/CONST/ G,P1,P12,RAD,HCONVR, SCONVR, DX,DY,DX2,DY2,1,0HEG. 
AM,N,M1,H2,NMl jNH2, IET, IDRY, IETP],CNTRANG,HEACT, DLEVEL 
COMMON/PRINTC/IP1, IP2, INC, JP1,JP2,DMULT, HMULT, ZMULT, RKNULT 
COMMON/CONS2/IGRID, IREF, ITAX, ITHAX, ID IEE 
c 


Crrkakkrkk COMPUTE THE WAVE NUMBERS USING THE PADE APPROXIMATION 


Chkkkkkkh (SEE THE COASTAL ENGINEERING NOTEBOOK) 
c 

DO 1 I=IWET,H 

DO 1 J=1,N 

DD=D(1,J) 
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DUMI=(OMEGAK2) ADD/G 
DUM2=DUM1+1.0/(1.0+0. 65224DUM1 +0. 46224 
K(DUK1AX2) +0. 08644 ( DUM A44)+0.06794 (DUM AAS) ) 
DUMS=TASQRT (GADD/DUM2 ) 
RKNUM=P 12/DUM3 
RKA( I,J) =RKNUM 
GRDK(1,J)=RKNUM 
C 
Carkkkakk COMPUTE THE HEIGHTS AND ANGLES USING SNELL’S LAW 
Cikarkakk COMPUTE SIGMA AND CCG 
C 
DUM1=RKNUMADD 
DUM2=2. OADUM] 
DUM3=TANH (DUML ) 
DUM4=S INH (DUM2) 
SINE=S IN (THETAO-CNTRANG) ADUM3 
ZANG=P I-ASIN(S INE) +CNTRANG 
2(1,J)=ZANG 
SI(I,J)=SIN(ZANG) 
CO(I,J)=COS(ZANG) 
DUHS=SQRT (COS ( THETAO-CNTRANG) /COS ( ZANG-CNTRANG) ) 
H(T,J)=HHADUMSASQRT(0.9/(0.5%(1.0+DUN2/DUM4) 
AKDUM3) ) 
COG(T,J)=0.54(1.0+(DUM2/DUM4) ) \ OMEGAR2) / 
A(RKNUMAK2) 
1 CONTINUE 
C 
Cakkrkakk CHECK FOR WAVE BREAKING 
C 
DO 600 I=IWET,M 
DO 604 J=1,N 
DUAS=H(1, J} 
HBRK=0.78AD(1,J) 
TE(DUMS.GE.HBRK) DUHS=HBRK 
C 
CAARARAARCONVERT HEIGHTS TO AMPLITUDE FUNCTION (HAG/2x5 IGNA) 
C 
H(1,J)=DUMS/HEACT 
604 CONTINUE 
600 CONTINUE 
RETURN 
END 
CAARRAARARAKRARARAARARAARARAARARRARRAARARRARRARRARARARARARRKARARKARAREARK 
C 
SUBROUTINE ANGLE( ISTART, IEND, ITMAX,WT,ALFHA, SCNV) 
C 
CAAKAAKAKKARKARAAKARAAARKARRARAARARAAARARARARARRARRARRARRARAARARARARAR AK 
C 
C THIS SUBROUTINE ITERATES TO SOLVE FOR THE WAVE ANGLES 
C ALONG A GIVEN ROW 
C 
CARRAKAAARAKARKARARARARARARARARARRKAARARRARAARAARAARARARARARRARARRRAR ARK 
C 
PARAMETER( [Q=95,, JQ=95) 
C 
COMMON/ANGLES/7Z(1Q, JQ) ,S1(10,JQ) .CO{ 10, JQ) 
COMMON/DEPTHS/D( 1a, JQ) 
COMMON/WAVNUM/RKA( 10, JQ) ,GRDK( 10, JQ) 
COMMON/CONST/ G,PI,P12,RAD,HCONVR, SCONVR,DX,DY,DX2,DY2,1,OMEG, 
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kM,N,M1,M2,NM1,NM2, IWET, IDRY, IWETP1 ,CNTRANG,HEACT, DLEVEL 
COMMON/CONS2/IGRID, IREE, ITAMX, ITHAX, IDIFE 
COMMON/MUSC/XMUC (IQ), XMUS( TQ), YAUC( JQ) , YMUS (JQ) 
DIMENSION DUM1( 10, JQ) ,X(JQ) 
C 
CrkkAAAAK SOLVE THE [IEFERENCED FORM OF THE IRROTATIONALITY OF THE 
Ckkkkkkkk GRADIENT OF THE WAVE PHASE FUNCTION EQUATION 
C 
DO 1 I=IEND, ISTART,-1 
IM=I-1 
LO 4 IT=1, ITHAX 
DO 2 J=2,NM1 
DUMG=ALPHAA(GRDK( 1, J+1)&S1(1,J+1)) 
k +(1.0-2. OAALPHA) A(GRDK(1,J)4S1(1,J)) 
k +ALFHAA(GRDK (1,J-1)AS1(T,J-1)) 
DUM2=WTA(GRDK(I-1,J+1)ACO( I-1,J+1)-GRDK(I-1,J-1)A 
RCO(I-1,J-1))/(DY2AYHUC (I) ) 
k +(1.0-WT) ACGRDK (1, J+1)4CO( 1, J+1)-GROK(1,J-1)% 
RCO(I, J-1))/(DY2RYMUC (I) ) 
TIUMS=DUMG- (DXAXMUS ( T-1) ) ADUN2 
DUM4=DUM3/GRDK< 1-1, J) 
DUM1(1,J)=DUM4 
2 CONTINUE 
SUK=0.0 
BO 22 J=2,NHl 
DUM4=DUM1(1, J) 
DWHS=Z( 1-1, J) 
DUK3=S IN(DUMS ) 


C 
CkkkkAKKA UPDATE THE SINES 
C 
EPSZ=0.9 
DUM4=EPSZADUMA+( 1. 0-EPSZ) ADUMS 
C 
Cakkkkhkk LIMIT OBLIQUE WAVE ANGLES 
C 


IF(DUM4.GE.0.997) DUM4=0.997 
TE(DUM4 .LE.-0.997) DUM4=-0.997 
DUM6=P 1-AS IN(DUMA4) 
SUM=SUM+ABS ( DUM6-Z( IK, J) ) 
2(1M,J)=DUM6 
SI(IM, J)=DUM4 
CO‘ 1M, J)=COS(DUMG) 
22 CONTINUE 

G 

ChakARAAA SET LATERAL BOUNDARY CONDITIONS FOR WAVE ANGLES ,SINES,AND 

Carkakrkk COSINES 


C 
Z(H, 1)=Z( 1M, 2) 
CO(IM,1)=CO( IK, 2) 
SIC IM,1)=S1( IM, 2) 
Z(TK,N)=Z( IM, NM1) 
CO(1M,N)=CO¢ IM, NM1) 
SICIM,N)=S1( IM, NHL) 

C 

ChkkkAAAK CHECK FOR ANGLE CONVERGENCE 

C 


TE(SUM.LT. (NM2ASCNV) )G0 TO 1 
4 CONTINUE 
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IE(IREE.EQ.1)G0 TO 1 ANGLE 
WRITE(6,643) IM ANGLE 
643 FORMAT(1X, ‘RELAXATION FOR WAVE ANGLES FAILED ON ROW =’, 14) ANGLE 
1 CONTINUE ANGLE 
99 RETURN ANGLE 
END ANGLE 
CRARAAKARAARARAARARARARARARRARRARAAAARARARRARRRARRARRARARAAARAAAARAAKARK DEPTH 
C DEPTH 
SUBROUTINE DEPTH DEPTH 
C DEPTH 
CAAAAKAKAAAAAAAAKARAAAARAAAARARARARRARAAARARRARARAAARAARRARAKARARKARRAKK DEPTH 
C DEPTH 
C THIS SUBROUTINE GENERATES THE WATER DEPTH MATRIX DEPTH 
C DEPTH 
CAKAAAAAARAAKARAARARAARAKAARARARRARARARARRARARRAARRARARRKRARAARAAAAARAAK DEPTH 
C PARAM 
PARAMETER ( 1Q=95, JQ=95) PARAM 
C PARAM 
COMMON/DEPTHS/D( IQ, JQ) DEPTH 
COMMON/CONST/ G,P1,P12,RAD,HCONVR, SCONVR,DX,DY,DX2,0Y2,1,0MEG, DEPTH 
AM,N,M1,M2,NM1,NM2, IWET, IDRY, IWETP1,CNTRANG,HEACT, DLEVEL DEPTH 
COMMON/CONS2/IGRID, IREF, ITAMX, ITHMX, IDIFE DEPTH 
COMMON/MUSC/XMUC (IQ) ,XMUS( 1), YMUC (JQ), YHUS(JQ) DEPTH 
COMMON/COOR/XX (10), YY(JQ) DEPTH 
C DEPTH 
CAkAKAKKK READ IN VARIABLE BATHYMETRY FROM FILE CODE 38. FIRST READ DEPTH 
CAkRAAAAK IN THE ’ALONGSHORE’ ROW CLOSEST TO SHORE, THEN PROCEED DEPTH 
CrkRAKKAK ROW BY ROW IN THE OFFSHORE DIRECTION DEPTH 
C DEPTH 
DO 14 I=1,4 DEPTH 
READ(8,90)(D(1,J),J=1,N) DEPTH 
14 CONTINUE DEPTH 
30 EORMAT(10F8.2) [TEPTH 
99 RETURN DEPTH 
END DEPTH 
CARARKARAARARARARARARARARARARRRARRARARARARARARARKKAKKARKKARKARRKARARARAAR POUT 
C POUT 
SUBROUTINE POUT(II1, 112, IYAL, JSTART, JEND, ITITLE, DUM1,FACT) POUT 
C POUT 
CARAKAAAARARKAAAKARARARARARARRARARARARARKARAAARARARARARARRARRRRRARARRARKA POUT 
Ct POUT 
C THIS SUBROUTINE PRINTS OUT SELECTED VALUES OF A PARTICULAR ARRAY POUT 
C ,DUM1, WHICH ARE SCALED BY THE VALUE OF ‘FACT’ POUT 
c POUT 
CARRAARRAARARARAAARARRARARARARARRRARARARRARAARARAARARARRRARRARRRARRARARR POUT 
C PARAM 
PARAMETER( IQ=95, JQ=95) PARAM 
C PARAM 
COMMON/CONST/ G,PI,PI2,RAD,HCONVR, SCONVR,DX,DY,DX2,D¥2,1,QMEG, POUT 
AM,N,H1,H2,NM1,NM2, IWET, IDRY, IWETP1 CNTRANG,HEACT , DLEVEL POUT 
COMMON/CONS2/ IGRID, IREF, ITAMX, ITHMX, IDIFE POUT 
INTEGER IX(JO+31), ITITLE(3) POUT 
DIMENSION DUM1( 18, JQ) POUT 
C POUT 
NC=31 POUT 
WRITE(6, 100) ITITLE, FACT POUT 
100 FORMAT(///,3A10,5X, ‘(MULTIPLIED BY ’,F6.0,’)’) POUT 
J1=ISTART POUT 
J2=ISTART+NC-1 POUT 
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1 IF(J2.GI.JEND) J2=JEND 
WRITE(6,102)(J,J=J1,J2) 
102 FORMAT(/,3X,’I/J:‘,3114) 
WRITE(6,103) 
103 FORMAT(1X, 


k 1 eee ee -  - - -  - -  - e - - - - -- - = = - = 


k 1 en nee enon saa owen My 


DO 2 I=L, 112, IVAL 
DO 3 J=J1,J2 
RND=0.9 
IF(DUM1(1,J).LT.0.0)RND=-0.5 
3 IX(J)=INT( EACTADUM] (1,J)+RND) 
2 WRITE(6,104) 1, (1X(J),J=J1,J2) 
104 FORMAT(1X,13,2X,‘:’,3114) 
J1=J1+NC 
J2=J2+NC 
IF(J1.LE.JEND)GO TO 1 
WRITE(6,101) 
101 FORMAT(//) 
RETURN 
END 


SUBROUTINE GRID 


THIS SUBROUTINE READS THE VARIABLE GRID MAPPING INFORMATION FROM 


AAKARKKARAKARKARRAAKAARARARRARRARARARRRRARRARARRARRARRRARRERRRARRRRRRKR 


FILE CODE 3 (X REGION FIRST, THEN Y) AND GENERATES THE VARIABLE 


GRID SCALE FACTORS ANI! DISTANCES TQ THE GRID CENTERS 


WHEN GENERATING A VARIABLY SIZED GRID IT IS RECCOMENDED THAT THE 
MAPPING BE DONE IN MAP INCEHES FROM SOME BATHYMETRIC CHART OF A 
KNOWN SCALE. THE ‘DX’ AND ‘DY’ USED AS INPUT INTO THE MODEL ARE 
THEN SIMPLY THE NUMBER OF FEET OR METERS PER MAP INCH. 


IF MAP 


INCHES ARE NOT USED IN THE MAPPING, THE VALUES OF ‘XSCALE’ AND 


‘YSCALE’ CAN BE SET 10 SOMETHING OTHER THAN 1.0 (IN CONJUNCTION 


WITH THE CHOICE OF ‘DX’ AND ‘DY’) TO CONVERT FROM THE MAPPELI 


SPACE INTO REAL SPACE 


KARKARAKARKARAARRAARARRRARARRARARARARARRARARRARARARARARARAARARARARARARKA 


ARDEFINITIONS OF VARIABLE GRID MAPPING INPUT REQUIRED BY THE MODEL 


NREG(X,Y) - IS THE NUMBER OF MAPPING REGIONS IN THE X OR Y 


DIRECTION 


NA(X,Y) - IS THE NUMBER OF GRID CELLS IN A PARTICULAR MAPPING 
REGION IN THE X OR Y DIRECTION 
A(X, Y)pB(X,%) ,C(X,¥) - ARE THE MAPPING COEFFICIENTS FOR A 


PARTICULAR REGION IN THE X OR Y DIRECTION 


PARAMETER ( IQ=95, JQ=95) 


KAKARKAARKAARARARRARRAARAARARARRARARRARRARARRARARRARAARRARRARRRRRARRRRR 


COMMON/CONST/ G,PI,P12,RAD,HCONVR,SCONVR,DX,DY,DX2,DY2,1,OHEG, 


AM,N,HL,M2,NML,NH2, IWET, IDRY, IWETP1,CNTRANG,HEACT, DLEVEL 


COMMON/CONS2/IGRID, IREF , ITAMX, ITHAX, IDIEE 


COMMON/HUSC/XMUC( TQ), XMUS( TQ), YMUC( JQ). YHUS( JQ) 
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C 


Caarrkarkk READ MAPPING INFORMATION 


C 


C 
CAkkARKAA WRITE THE VARIABLE GRID INFORMATION ON FILE CODE 6 


C 


COMMON/COOR/XX( TQ), YY (JQ) 


DIMENSION AX(50),BX (50) ,CX (SO) .NAX(S0) 
DIMENSION AY(50),BY(50) ,CY(50) NAY (50) 


CHARACTERA] SIR 


XSCALE=1.0 
YSCALE=1.0 


I=] 
I=] 
2] CONTINUE 


READ(3,31,END=999) STR,NST,NEND,Al,B1,C1 


31 FORMAT(AL,7X,218,3016.11) 
NEND=NEND-1 
IE(STR.EQ.’Y’) GO TO 20 
AX(1)=Al 
BX(1)=1 
CX(1)=C1 
DO 201 LLC=NST,NEND 
LL=LLC+1 
ALPHC=LL-0.5 
ALPHS=LL 


XMUS (LL-1)=(B1AC1A(ALPHSAR(C1-1.0)) )AXSCALE 
XMUC(LL-1)=(BIAC1A(ALPHCAR(C1-1.0)) )AXSCALE 


XX(LL-1)=(AL+B1A(ALPHCAACI ) )AXSCALEALX 


201 CONTINUE 
NAX(1)=NENI+1-NST 
NXREG=1 
I=I+1 
G0 TO 21 

20 CONTINUE 
AY(J)=Al 
BY(J)=B1 
CY¥(J)=C1 
10 203 LLC=NST,NEND 
LL=LLC+1 
ALPHC=LL-0.5 
ALPHS=LL 


YMUS (LL-1)=(BIACLA(ALPHSAA(C1-1.0))) AYSCALE 
YHUC (LL-1)=(B1AC1&(ALPHCAR(C1-1.0)) }&YSCALE 


YY(LL-1)=(A1+B1A(ALPHCAACL ) ) AYSCALEALY 


203 CONTINUE 
NAY (J) =NEND+1-NST 


999 CONTINUE 


WRITE(6,9) 


9 EORMAT(////,20X, VARIABLE WAVE GRID INFORMATION’ ,///) 


WRITE(G,23) XSCALE, YSCALE 


23 EORMAT(///, EXPANSION COEFFICIENT SCALE FACTORS’./, 


k’XSCALE=’ ,F10.3,0X, ‘ YSCALE=’ ,F10.3,//) 


WRITE(6,9) 


5 FORMAT(1X,°X EXPANSION COEFS(A,B,C) AND GRIDS PER REGION’ ,//) 
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WRITE(6,4)¢1,AX(1) BACT) ,CX(1) ,MAX(1), I=1,NXREG) GRID 


4 FORMAT(110,3E20.7, 110) GRID 
WRITE(6,6) GRID 

6 FORMAT(///,1X,’Y EXPANSION COEFS(A,B,C) AND GRIDS PER REGION’,//) GRID 
WRITE(6,4)(1,AY(1),BY(1) CY(1) MAY (1), I=1,NYREG) GRID 
WRITE(6,7) GRID 

7 FORMAT(///,1X,’X SCALE FACTORS (XMUC,XMUS)’,//) GRID 
WRITE(6,2)(1,XMUC( I) ,XMUS(I), I=1,H) GRID 
WRITE(6,8) GRID 

® EORMAT(///,1X,’Y SCALE FACTORS (YMUC,YHUS)’,//’) GRID 
WRITE(6,2)(1,YMUC(1), YMUS( 1), I=1,N) GRID 

2 FORMAT(S(15,2F10.9)) GRID 
WRITE(6,10) GRID 

10 FORMAT(///,1X,‘X CENTER DISTANCES’ ,//} GRID 
WRITE(6,11)(1,XX(1),1=1,4) GRID 

11 FORMAT(8(14,F12.3)) GRID 
WRITE(6,12) GRID 

12 FORMAT(///,1X,’Y CENTER DISTANCES’ ,//) Gk ID 
WRITE(G,11)(1,YY¥(J),J=1,N) GRID 
RETURN GRID 

END GRID 
CARAKARAAARARARAAARAAAKARAR RARAAAARARAARARRARARARARARARRARRARRRARRRARAKK ONBC 
C ONBC 
SUBROUTINE ONBC(DUM1) ONBC 

C ONBC 
CARARAAAKARRAAKARARARAAARARRARARAARARRRRARARARARAAAAKAARAARAARRARARKARKK ONBC 
C ONBC 
C THIS SUBROUTINE SETS THE ONSHORE BOUNDARY CONDITION . ONBC 
C THIS BOUNDARY CONDITION IS NOT USED IN THE COMPUTATIONAL SCHEME ONBC 
C SINCE IT PROCEEDS FROM OFFSHORE TOWARDS ONSHORE. ONBC 
c ONEC 
CARRARAARAAAAAAAARARARAAAARAAAARARRRARRARARRARARARAAARAARKAARARAARAARARR ONBC 
C PARAM 
PARAMETER ( 10=95, JQ=95) PARAM 
C PARAM 
COMMON/CONST/ G,PI,P12,RAD,HCONVR,SCONVR DX ,DY,DX2,DY2,1,OMEG, ONBC 

AHN ML M2,NM1.NM2, IWET, DRY, IWETP ,CNTRANG HFACT, DLEVEL ONBC 
COMMON/CONS2/ IGRID, IREF, ITAMX, ITHAX, IDIFE ONBC 
DIMENSION DUM1( 10, JQ) ONEC 

0 1 I=1, IDRY ONBC 

DO 1 J=1.N ONBC 
TWH1(1,J)=DUM] (IWET, J) ONBC 

1 CONTINUE NEC 
RETURN ONEC 

END ONBC 
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APPENDIX G: SAMPLE FILES--FIELD RESEARCH FACILITY PIER, 
DUCK, NORTH CAROLINA, WAVE PROPAGATION EXAMPLE 
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FIELD RESEARCH FACILITY (DUCK , NC) EXAMPLE 


re) 12.0 50 24.0 9.800 0.0 2 0.90 
2.00 12.00 20.00 
1.50 8.00 -35.00 

oar!) Dale b Der} 
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Figure G3. DUCKDAT 


ID MODS 
D PARAN. 3 
PARAMETER (1G=751 JQ=50) 
D DEPTH. 23; DEPTH. 26 
DO 14 J=1,N 
READ (8,60) (D(1,J).1=15M) 
DO 315 I=1,8 
D(1,J)=-D(1,J) 
315 CONTINUE 
60 FORMAT (SX/(15F8.2)) 
14 CONTINUE 
D MAIN, 203 
DMULT=10.0 


Figure G4. DUCKUPD 


ID MODS 
D PARAM. 3 
PARAMETER (1@=75, JQ=30) 
D DEPTH. 23, DEPTH. 26 
DO 14 JF1.N 
READ (8460) (B(1, J) 5 f=154) 
DO 315 I=1,4 
D1, =-D(1 J) 
315 CONTINUE 
60 FORMAT (SX/ (158. 2) ) 
14 CONTINUE 
D MAIN. 203 
DNULT=10.0 
D MAIN.2 
PROGRAM RCPWAVE (INPUT, OUTPUT, TAPE7=TAPE7» TAPES=TAPE6, TAPEB=TAPES 


%*, TAPES=TAPE3) 
Figure G5. DUCKUP2 
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REGIONAL COASTAL PROCESSES WAVE TRANSFORMATION MODEL 
(RCPWAVE) 


FIELD RESEARCH FACILITY (DUCK » NC) EXAMPLE 


MODEL INPUT: 


UNIFORMLY SIZED, RECTILINEAR GRID MESH 


X - DIRECTION (ON-OFFSHORE) 73 CELLS, EACH 12.00 IN LENGTH 
Y - DIRECTION (LONGSHORE) 30 CELLS, EACH 24.00 IN LENGTH 


THE ACCELERATION OF GRAVITY IS 9.80. THESE UNITS DETERMINE THE UNITS OF ALL INPUT AND OUTPUT VARIABLES 


20.000 
- 35.000 


THE DEEP WATER WAVE PARAMETERS FOR CASE 1 ARE: HEIGHT= 2.000 PERIOD= 12.000 ANGLE 
THE TIEEF WATER WAVE PARAMETERS FOR CASE 2 ARE: HETGHT= 1.500 PERIOD= 9.000 ANGLE 


THE OFFSHORE CONTCURS MAKE AN ANGLE OF = .00 DEGREES WITH THE Y AXIS 
THE BATHYMETRY MATRIX 1S VARIABLE IN BOTH HORIZONTAL DIRECTIONS AND WAS READ FROM FILE CODE 8 


A WATER LEVEL CHANGE OF  .90 WAS ADDED TO THE ENTIRE BATHYMETRY MATRIX 


Figure G7. RCPPRNT 
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THE DEEP WATER WAVE PARAMETERS FOR CASE 1 ARE 


WAVE CONDITION 
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Figure G8. 
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WAVE NUMBER (MULTIPLIED BY 1000.) 
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HOMER SPIT , ALASKA EXAMPLE 
96 416.7 83 833.3 32.200 0:05 oa 0.00 
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130 1 28 0% 


0 
Figure H3. HOMRDAT 
*ID MODS 
*D PARAM.3 
PARAMETER( 10=96 , JQ=83) 
a] DEPTH.26 
DO 10 I=1,M 
DO 10 J=1,N 
D(1,J)=D(1 ,J)*6.0 
10 CONTINUE 
Figure H4. HOMRUPD 
*1D MODS 
*D PARAM.3 
PARAMETER(10=96, JQ=83) 
x1 DEPTH. 26 
DO 10 1=1,M 
DO 10 J=1,N 
D(1,J)=D(1 ,J)*6.0 
10 CONTINUE 
*D MAIN.2 


PROGRAM RCPHAVE( INPUT , OUTPUT , TAPE7=TAPE7 , TAPE6=TAPE6 , TAPES=TAPES 
+, TAPE3=TAPE3) 


Figure H5. HOMRUP2 
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Figure H6. HOMRDEP 
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REGIONAL COASTAL PROCESSES WAVE TRANSFORMATION MODEL 
(ROPWAVE) 


HOMER SPIT » ALASKA EXAMPLE 


MODEL INPUT: 


UNIFORMLY SIZED, RECTILINEAR GRID MESH 


X - DIRECTION (ON-GFFSHORE) 96 CELLS, EACH 416.70 IN LENGTH 
Y - DIRECTION (LONGSHORE) 83 CELLS; EACH 833.30 IN LENGTH 


THE ACCELERATION OF GRAVITY 1S 32.20. THESE UNITS DETERMINE THE UNITS OF ALL INPUT AND OUTPUT VARIABLES 
THE DEEP WATER WAVE PARAMETERS FOR CASE 1 ARE: HEIGHT= 8.000 PERTOD= 10.000 ANGLE= -45.000 
JHE OFFSHORE CONTOURS MAKE AN ANGLE OF  § .00 DEGREES WITH THE Y AXIS 

THE BATHYMETRY MATRIX IS VARIABLE IN BOTH HORIZONTAL DIRECTIONS AND WAS READ FROM FILE CODE 8 


A WATER LEVEL CHANGE OF  .00 WAS ADDED TO THE ENTIRE BATHYMETRY MATRIX 
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WAVE CONDITION 
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-14 -13 -13 -14 -15 -17 -19 -21 -21 -21 -23 -25 -27 -27 -27 -29 -30 -30 -29 -28 -27 -29 -29 -27 -23 
33: -14 -14 -14 -14 -15 -17 -19 -21 -21 -21 -24 -25 -27 -2B -29 -30 -31 -30 -29 -28 -28 -29 -29 -26 -23 
34: -14 -14 -14 -15 -16 -18 -20 -21 -21 -22 -24 -26 -28 -30 -31 -32 -32 -31 -29 -28 -28 -29 -28 -26 -23 
: 14-14 -14 -15 -17 -18 -21 -22 -22 -24 -26 -28 -30 -32 -32 -33 -32 -31 -29 -28 -29 -30 -29 -26 -23 
36 : -13 -14 -15 -16 -17 -19 -21 -23 -24 -26 -27 -29 -31 -32 -33 -33 -32 -31 -27 -27 -30 -30 -27 -26 -22 
37: -14 -14 -16 -18 -19 -20 -22 -24 -25 -27 -28 -30 -31 -33 -33 -33 -32 -30 -29 -28 -29 -30 -28 -25 -21 
38: -14-15 -17 -19 -20 -21 -24 -25 -27 -28 -29 -3i -32 -34 -33 -32 -31 -30 -28 -28 -29 -27 -27 -24 -22 
39 : -15 -16 -18 -20 -21 -23 -25 -27 -28 -29 -30 -31 -32 -33 -33 -32 -31 -30 -28 -28 -29 -29 -27 -24 -22 
= -16 -16 -18 -20 -22 -24 -26 -28 -29 -30 -30 -31 -33 -33 -32 -32 -31 -30 -28 -28 -27 -28 -26 -24 -22 
Al : -17 -17 -19 -21 -23 -25 -27 -30 -31 -31 -31 -31 -33 -33 -32 -31 -30 -29 -28 -27 -28 -27 -25 -23 -21 

2 $ -17 -18 -21 -23 -24 -26 -29 -31 -32 -31 -30 -30 -31 -32 -31 -30 -30 -29 -27 -27 -27 -27 -25 -23 -20 
43: -18 -20 -22 -24 -25 -27 -30 -32 -32 -30 -29 -29 -30 -31 -30 -30 -30 -29 -27 -26 -26 -26 -25 -23 -20 
© -19 -21 -23 -25 -26 -28 -31 -32 -31 -30 -28 -28 -29 -30 -30 -30 -29 -28 -27 -26 -26 -25 -24 -21 -19 
45: -21 -23 -25 -26 -27 -29 -31 -32 -31 -29 -28 -28 -30 -31 -30 -30 -29 -28 -26 -26 -25 -24 -22 -20 -19 
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APPENDIX I: INTPRCP PROGRAM LISTING 


ooeos 


sees 


PotD 


sCOMDECK PARAM 


C 
PARAMETER (1Q=95, JQ=95) 
PARAMETER (182=95, JQ2=95) 
PARAMETER (KQ=18% JQ) 
C 
sDECK HAIN 
PROGRAM GRID(TAPE1=TAPE1, TAPE2=TAPE2, TAPEI=TAPES, TAPE4=TAPE4, 
LTAPEG=TAPE6, TAPE7=TAPE7) 
sCALL PARAM 
Ceeeesececceeercesscescecessenceceeeretececcesteeceecessces 
C 
C DOCUMENTATION 
C 
Coeeseeeaegeserrensessenssgaseseeesesesneseteses sees seness 
C 


C THIS PROGRAM PROVIDES A METHOD FOR INTERPOLATING BATHYMETRY 
C FROM ONE GRID TO ANOTHER. THE NEW GRID ORIGIN CAN BE 
C TRANSLATED AND/OR ROTATED RELATIVE TO THE OLD GRID ORIGIN 


C 
CeoscesereereeceesecreceeeereeeesereoeeetererereeseseeEeTEe 
C 

C INPUT DATA (TAPE7) 

C 


Cesescegegeeeresageccecrccseacegeserseseeeseceeeeseeeeeeesse 
C 


op] 


READ (7,500) HN» M2.N2,ANG, XSHIFT, YSHIFT 
FORMAT (415, 3F 10.5) 


$ 


M= NUMBER OF GRID CELLS IN THE X-DIRECTION ON THE OLD GRID 
N= NUMBER OF GRID CELLS IN THE Y-DIRECTION ON THE OLD GRID 
M2= NUMBER OF GRID CELLS IN THE X-DIRECTION ON THE NEW GRID 
N2= NUMBER OF GRID CELLS IN THE Y-DIRECTION ON THE NEW GRID 


ANG=THE ANGLE OF ROTATION MEASURED 
FROA THE OLD GRID TO THE NEW GRID 
+COUNTER-CLOCK WISE 
-CLOCKWISE 


XSHIFT=THE X OFFSET FROM THE OLD GRID TO THE NEW GRID 
+ NEW GRID ORIGIN IS BELOW THE OLD ORIGIN 
- NEW GRID ORIGIN IS ABOVE THE OLD ORIGIN 
YSHIFT=THE Y OFFSET FORM THE OLD GRID TO THE NEW GRID 
+ NEW GRID ORIGIN 1S TO THE RIGHT OF THE OLD ORIGIN 
- NEW GRID ORIGIN IS TO THE LEFT OF THE OLD ORIGIN 


sees 


READ(7,505) GRDTYP1,GRDTYP2 
505 FORMAT (215) 


GRDTYP1=0LD GRID TYPE 
0 = CONSTANT SIZED, RECTILINEAR 


ANINIONININANANANANANYNNONIANYNYANNAANMDNNYMNONAAIAGIHA 


T2 


Onrororwnre For WN 


32 


EaSZFSESSYF HEH 


= VARTABLY SIZED, RECTILINEAR 


GROTYP2=NEW GRID TYPE 


C 

C 

C 

C 0 = CONSTANT SIZED, RECTILINEAR 
C 1 = VARIABLY SIZED, RECTILINEAR 
C 

Cease 

C 

C IF: GRDTYP1=0 THEN: 

Cc 

C READ (7,501) DX,DY 

C S01 FORMAT (2F 10.5) 

C 

Coens 

C 

C IF: GRDTYP2=0 THEN: 

C 


C READ(7,501) DX2,DY2 
C S01 FORMAT (2F 10.5) 


C 

Ceeee 

C 

Cc IF: GRDTYP1=0 & GRDTYP2=0 THEN: 

C 

C READ(7,501) DX,DY 

C READ(7,501) DX2,DY2 

C 501 FORMAT (2F10.5) 

C 

C DX=CELL LENGTH IN THE X DIRECTION (OLD GRID) IN MAP INCHES OR CH 
C DY=CELL LENGTH IN THE Y DIRECTION (OLD GRID) IN MAP INCHES OR CH 
C DX2=CELL LENGTH IN THE X DIRECTION (NEW GRID) IN MAP INCHES OR CH 
Cc DY2=CELL LENGTH IN THE Y DIRECTION (NEW GRID) IN MAP INCHES OR CH 
C 

ATTETT Tie tittitiittititetiits tite i irri ti itetiitiisi itis 

C 

C INPUT DATA (TAPES) 

C -OLD GRID- 

Cc 

TITS TIT Ti titers tri t tsetse tt ititititrtir rier tier ere 

C 

C esNOTEse 

C ONLY USE TAPES IF GRDTYP1=1 

C 

COMER SREROAHETSSHERORAEATES ARS KETHETASEHTTSHERE STE RES ETORTS 


IF GRDTYP1=1 THEN: 


READ (3,31, END=999) STR»NST,NEND,AI,B1,C1 


STR=X OR Y STRETCHING DIRECTION 

NST=STARTING CELL NUMBER FOR THE REGION 
NEND-NST=NUMBER OF CELLS IN THE REGION 
A1,Bi,C1=STRETCHING COEFFICIENTS FOR THE REGION 


ooonrnaninniwnnn 


13 


MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
RAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
HAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 


SBGEFARSARSAESSSIRFHLASE SSS 


C 

C 31 FORMAT(AL, 7X, 218, 3616.11) 

C 
CRMREKSSRSHTEREKESESKESEEKEKEEECEREDT SRT ESERESESEREE STEER EE 
C 

C INPUT DATA (TAPE4) 

C -NEW GRID- 

C 
COOGERGEEETORTETTNGTEDNGDSONTNED EST GETR SDDS TFSSS 9999S OTE SS 
C 

C soNOTEse 

C ONLY USE TAPE4 IF GRDTYP2=1 

C 


CeSeegaageeeKee see ESE HTS SSSS KE SLSCSESSSSTLESRESSSS HS ESASSSES 


C IF GRDTYP2=1 THEN: 
READ (4,31,END=778) STR»NST,NEND,A1,B1,Ci 


C 
C 
C 
C STR=X OR Y STRETCHING DIRECTION 

C NST=STARTING CELL NUMBER FOR THE REGION 

C NEND-NST=NUMBER OF CELLS IN THE REGION 

C Al,B1,C1=STRETCHING COEFFICIENTS FOR THE REGION 
C 

C 

C 

C 

C 

Cc 


31 FORMAT {AL, 7X, 218, 3616. 11) 


CHCKHKHLH KAAS SHAS SKK HSE HSH RHSSTESS SLES ER ESE SERRE 


C 
C INPUT DATA (TAPE1) 
C 
C 


SAGCHRHRAHHORESS ASS LSVSHTAKIH SHAK SSHHSS RSS HHLAKGHSSSSHS SSS OS 


C 

C DO 10 I=1,4 

C READ(1,11) (DD(1,J),J=1.N) 
C 10 CONTINUE 

C 11 FORMAT (1OF8. 2) 
C 

C 

C 

C 


DD(1,J)=DEPTH AT EACH CELL OF THE OLD GRID 


CHRMAROTHAAERERAORSESLSKRRERSRERSKO TET SRE TEESE SSSEST SETAE THs Es 
COMMON/CONST/M;>N,M2.N2, DX, DY 
COMMON/COOR/XX (18) YY (JQ) 
COMMON /COGR2/XX2 (182) » Y¥2(J@2) 
COMMON/ COORI/XOMRTN (1, JQ), YOWRTN (18; JQ) 
COMMON/DEPTHS/D(1@2,J@2) ,DD(1Q, JQ) 
CHARACTER#1 STR 
READ (7590) MyNsM2,N2,ANG) XSHIFT> YSHIFT 
ANG=-1. 0# (ANG) 
READ(7,505) GRDTYP1, GRDTYP2 

305 FORMAT (215) 

500 FORMAT (415, 5F 10.5) 
ANG=ANG#3. 1415927/180.0 
SCALE=1.0 
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MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
PAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
BAIN 
MAIN 
MAIN 
MAIN 
MAIN 
HAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
RAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
RAIN 
HAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 


XSHIFT=XSHIFT*SCALE 

YSHIFT=YSHIFT#SCALE 

IF(GROTYP1.E@.0) GO TO 304 
21 CONTINUE 

READ (3,31, END=999) STR» NST, NEND»A1,B1,C1 
31 FORMAT‘A1, 7X, 218, 3616.11) 

NEND=NEND-1 

IF(STR.EQ.’Y’) GO TO 20 

DO 201 LLC=NST.NEND 


XX(LL-1)=(A1+B1¢ (ALPHC##C1) ) sSCALE 
201 CONTINUE 
60 TO 21 
20 CONTINUE 
DQ 203 LLC=NST,NEND 
LL=LLC+1 
ALPHC=LL-0.5 
ALPHS=LL 
YY(LL-1)=(A1+B19 (ALPHC##C1) ) #SCALE 
203 CONTINUE 
GO TO 21 
999 CONTINUE 
506 CONTINUE 
WRITE (619) 


9 FORMAT (////120X," VARIABLE WAVE GRID INFORMATION’ ///) 


WRITE (6123) SCALE 


23 FORMAT (///, EXPANSION COEFFICIENT SCALE FACTOR’, /, 


we SCALE=",F10.3,//) 
IF(GRDTYP1.E@.0) GO TQ 307 
G0 TO 513 
507 READ(7,501) DX,DY 
501 FORMAT ‘(2F10.5) 
DO S11 L=1,¥ 
XX(L)=(L-0. 5) #DXsSCALE 
S11 CONTINUE 
DO 512 L=i,N 
YY(L)=(L-0. 95) #DY#SCALE 
512 CONTINUE 
513 CONTINUE 
WRITE (6) 10) 
10 FORMAT (///,1X,7X CENTER DISTANCES’ »//) 
11 FORMAT‘(8(14,F12.3)) 
WRITE (6.12) 

12 FORMAT(///,1X,"Y CENTER DISTANCES’, //) 
WRITE (4,11) (J, YY (J). J=15N) 
IF(GRDTYP2.E8.0) GO TO 539 

28 CONTINUE 
READ (4, 31,END=998) STR» NST, NEND,A1,B1,C1 
NEND=NEND-1 
IF(STR.E@.'Y*) GO TO 27 
DO 271 LLC=NST.NEND 
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MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
HAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 


6oge% 


LL=LLC+1 

ALPHC=LL-0.5 

ALPHS=LL 

XX2(LL-1)=(A1+B1 (ALPHC#*C1) ) #SCALE 
271 CONTINUE 

GO TO 28 
27 ~CONTINUE 

DO 273 LLC=NST,NEND 

LL=LLC+1 


YY¥2(LL-1)=(A1+B1# (ALPHC##C1) ) #SCALE 
273 CONTINUE 
GO TO 28 
998 CONTINUE 
G0 TO 560 
539 CONTINUE 
READ(7,599) DX2,DY2 
J00 FORMAT (2F10.5) 
369 CONTINUE 
WRITE (4, 19) 
19 FORMAT (////,20X,°VARIABLE SECOND GRID INFORMATION’, ///) 
WRITE (6123) SCALE 
IF(GRDTYP2.E@.0) 60 TO 515 
60 TO 525 
S15 CONTINUE 
DO 520 L=1,42 
XX2(L)=(L-0.5) #DX2¢SCALE 
920 CONTINUE 
DO 524 L=1,N2 
YY2(L)=(L-0.5) #DY2*SCALE 
524 CONTINUE 
520 CONTINUE 
WRITE (6,40) 
AQ FORMAT (///,1X,'X2 CENTER DISTANCES’, //) 
WRITE (6.11) (1,XX2(1) » T=1,M2) 
WRITE (6,42) 
42 FORMAT (///,1X,°Y2 CENTER DISTANCES’, //) 
WRITE (6,11) (Js Y¥2(J) » J=1,N2) 
DO 26 JJ=1,N 
DO 246 II=1,h 
XOWRTN (11s JJ)=(XX (11) -XSHIFT) #COS (ANG) + (YSHIFT 
1-YY (JJ) ) #SIN (ANG) 
YOWRTN (11s JJ)=(XX (LL) -XSHIFT) #SIN (ANG) + (YY (JJ) 
1-YSHIFT) #COS (ANG) 
26 CONTINUE 
CALL INTERPO 


99 STOP 
END 
sDECK INTERPO 


CHSSOKHHEAHOATHHHLAARSHLRARTATAHESASORASHRSRS HSH SSS TASS AS RAVAGE TESS 


c+ 
SUBROUTINE INTERPO 
Ce 


16 


231 


wean BG BRASH 


PEE 


CARMA AKA A REAR AR ATA SHER AH ERA AHHK HRA AVAAR ARH HKO RAH AG SHREK HS HRS 


#CALL PARAM 


10 
{1 


579 


60 


pone 


COMMON /CONST/M,N,M2,N2,DX.DY 
COMMON/COOR/XX (18), YY (JQ) 

COMMON/ COOR2/XX2 (182). YY2(J@2) 
COMMON/COORS/XOWRTN (18; JQ) , YOWRTN (12, JQ) 
COMMON/DEPTHS/D(1@2, J@2) ,DD(1@, JQ) 
DIMENSION R(KQ)s IVAL (K@) » JVAL (KQ) , DIST (K@) 
INTEGER MD2,ND2 

DO 10 I=1." 

READ(1,11) (DD(I,J),J=1,N) 

CONTINUE 

FORMAT (10F8. 2) 

CALL POUT(1.M,1,1.N,°OLD DEPTHS (X.1)',DD,1.0,1@) JQ) 
DO i [=1,2 

DO 1 J=1,N2 

XD=XX2(1) 

YD=YY2 (J) 

K=0 

IF(1.£Q@.1.AND.J.EQ.1) G0 TO 401 

IF(J.NE.1) GO TO 350 

IiFP=1 1FR+30 

DO 599 II[=11FM, 1iFP 

IF(II.LE.0.0R.11.67.M) GO TO 599 

DO 5 JJ=1.N 

K=K+1 

DIST (K) =SRT ( (XOWRTN(TI, JJ)-XD) #¢2+(YOWRTN (IT, JJ)- 
1YD) ##2) 


CS=XOWRTN(II, JJ)-XD 
SN=YOMRTN (IT, JJ)-YD 
IF(CS.GE.0.0. AND. SN. GE.0.0)R(K)=1.0 
IF(CS.LT.0.0. AND. SN. GE.0.0)R(K)=2.0 
IF(CS.LT.0.0. AND. SN.LT.0.0)R(K)=3.0 
IF (CS. GE.0.0. AND. SN.LT.0.0)R(K)=4.0 
IF(CS.E@.0.0,AND.SN.E@.0.0) 60 TO 403 
IVAL(K)=I1 

JVAL (K)=JJ 

IF(K.GT.1) GO TO 565 

DHIN=DIST (1) 

T1=1VAL (1) 

Ji=IVAL (1) 

K1=1 

IF(DIST(K).GE.DMIN) 60 T0 5 
DHIN=DIST (K) 

kik 

CONTINUE 

CONTINUE 

T1=IVAL(K1) 

J1=JVAL(K1) 

Ri=R(K1) 

T1FM=11-15 

60 TO 552 

MD2=(H#1.0) 

ND2=(N#1.0) 


I7 


GOYPUSUBEESESESRGSSSSRUHHS HHS SRBYRBRRARVSSSGERSAGSES 0 


DO 602 [[=1,MD2 
DO 602 JJ=1,ND2 
K=K+1 
DIST (K) =SQRT ( (XOWRTN (II, JJ)-XD) ##2+ (YORRTN ( 
111, JJ)-YD) ##2) 
CS=XOWRTN (II, JJ)-XD 
SN=YORRTN (IT, JJ)-YD 
IF (CS. GE. 0.0. AND. SN. GE.0.0)R(K)=1.0 
IF(CS.LT.0.0. AND. SN. GE.0.0)R(K)=2.0 
IF(CS.LT.0.0. AND. SN.LT.0.0)R(K)=3.0 
IF (CS. GE.0.0. AND. SN.LT.0.0)R(K)=4.0 
IF(CS.E@.0.0.AND.SN.E@.0.0) 60 TO 403 
IVAL(K)=I1 
JVAL (K) =JJ 
IF(K.GT.1) GO TO 610 
DMIN=DIST (1) 
I1=I1 
Ji=JJ 
K1=1 
610 IF(DIST(K).GE.DMIN) GO TO 602 
DHIN=DIST (K) 
K1=K 
602 CONTINUE 
T1=IVAL(K1) 
J1=JVAL (KL) 
R1=R(K1) 
IFH=I1-15 
G0 TO 352 
T1=IVAL(K) 
J1=JVAL(K) 
T1FM=I1-15 
609 DiI, J)=DD(11, 31) 
60 TO 717 
613 SUM1=DMIN+DMIN2 
SUM2=SUM1 /DMIN+SUM1 / DHIN2 
SF1=SUM1 / (DHINSSUM2) 
SF2=SUM1/ (DMIN2*SUF2) 
D(1, J)=SF1DD (11, JL) +SF2*DD (12, J2) 
G0 TO 717 
614 SUM1=DAIN+DAIN2+DMINS 
SUM2=SUR1/DAIN+SUM1 /DMIN2+SUM1 /DAINS 
SF 1=SUM1/ (DMIN®SUM2) 
SF2=SUM1 / (DMIN2*SURZ) 
SFI=SUN1 / (DMINS*SUR2) 
D(1, J)=SF1sDD (11) Ji) +SF2*DD (12, J2)+SF3eDD (13, J3) 
60 TO 717 
530 TiM4=11-4 
TiP4=11+4 
JiM4=J1-4 
JiP4=J1i+4 
DO 551 I[=11K4,11P4 
IF(I1.LE.0.OR.11.G7.4) 60 TO S51 
DO 553 JJ=Jin4, JiP4 
IF(JJ.LE.0.0R.JJ.GT.N) G0 TO 353 


60: 


wa 


18 


78 


360 


352 


700 


701 


370 


K=K+1 

DIST (K)=SQRT ( (XOWRTN (II, JJ)-XD) ##2+ (YOWRTN (11, JJ)-YD) ##2) 
CS=XOWRTN (IT, JJ)-XD 
SN=YORRTN (TT, JJ)-YD 

IF (CS. GE.0.0. AND. SN. GE.0.0)R(K)=1.0 
IF(CS.LT.0.0. AND. SN. GE.0.0)R(K)=2.0 
IF (CS.LT.0.0. AND. SN.LT.0.0)R(K)=3.0 
IF(CS.GE.0.0. AND. SN.LT.0.0)R(K)=4.0 
IF(CS.EQ.0.0.AND.SN.E@.0.0) 60 TO 403 
IVAL(K)=IT 

JVAL (K)=JJ 

IF(K.GT.1) 60 TO 560 

DMIN=DIST (1) 

T1=IVAL (1) 

Ji=JVAL (1) 

K1=1 

IF(DIST(K).GE.OMIN) GO TO 553 
DNIN=DIST (K) 

K1=K 

CONTINUE 

CONTINUE 

T1=IVAL(K1) 

J1=JVAL (K1) 

RI=R(K1) 

CONT INUE 

DO 700 ILI=1,K 

R2=R(11T) 

IF(R2.E@.R1) GO TO 700 

K2=III 

I2=IVAL(ITI) 

J2=JVAL (IIT) 

DMIN2=DIST (III) 

G0 TO 701 

CONTINUE 

IF(R2.E@.R1) GO TO 609 

CONTINUE 

DO 570 L=I1I,K 

IF(L.E@.K1) 60 TO 570 
IF(R(L).EQ.R1) GO TO 570 
IF(DIST(L).GE.DMIN2) GO TO 570 
DMIN2=DIST(L) 

K2=L 

CONTINUE 

R2=R (K2) 

I2=IVAL (K2) 

J2=JVAL (K2) 

DO 703 JJJ=1,K 

RI=R (JIJ) 

IF(R3.EQ.R1.0R.R3.EQ.R2) GO TO 703 
K3=JJJ 

T3=IVAL(JJJ) 

Ja=JVAL (JIJ) 

DHINZ=DIST (JJJ) 

GQ TO 704 
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703 


704 


CONTINUE 

IF (R3.EQ.R2.0R.R3.EQ.R1) G0 TO 413 
CONTINUE 

DO 580 LL=JJJ+K 

IF(LL.E@.K1.OR.LL.E@.K2) 60 TO 380 
IF(R(LL).EQ.R1.0R.R(LL).EQ.R2) GO TO 580 
IF(DIST(LL).GE.DMIN3) GO TO 580 
DMIN3=DIST (LL) 

K3=LL 

CONTINUE 

T3=IVAL(K3) 

J3=JVAL (K3) 

RS=R (KS) 

DO 388 ILI=1,K 

R4=R (111) 

IF (R4.EQ.R3.OR.R4.E@.R2.0R.R4.E8.R1) 6&0 TO 588 


K4=I11 

T4=IVAL (IIT) 
J4=JVAL (111) 
DHIN4=DIST (III) 
60 TO 589 


588 CONTINUE 


IF (R4.EQ.R3.OR.R4.£0.R2.0R.R4.E@.R1) GO TO 414 


389 CONTINUE 


DO 590 LLL=III.K 

IF(LLL.E@.K1.0R.LLL.E@.K2) GO TO 590 

IF(LLL.E@.K3) G0 TO 590 

IF (DIST (LLL). GE.DMIN4) GO TO 590 

IF(R(LLL) .E@.R1.OR.R(LLL).E@.R2.0R.R(LLL) .EQ.R3) GO TO 570 
DMIN4=DIST (LLL) 

K4=LLL 


570 CONTINUE 


T4=IVAL (K4) 
J4=JVAL (K4) 

R4=R (K4) 

D1=DAIN 

D2=DMIN2 

DI=DMINS 

D4=DMIN4 
SUM1=D1+D2+D3+D4 
SUM2=SUN1 /D1+SUM1/D2+SUM1 /D3+SUN1 /D4 
SF1=SUM1/ (D1sSUM2) 
SF2=SUM1/ (D2*SUM2) 
SFI=SUM1/ (D3sSUMZ) 
SF4=SUM1/ (D4*SUM2) 


D(1,J)=SF1eDD (11, JL) +SF2eDD (12, J2)+SF3eDD (13, J3) +SF4eDD (14, J4) 


717 CONTINUE 


1 CONTINUE 
CALL POUT(1,M2.1,1,N2, NEW DEPTHS (X,1)’,D,1.0,1@2, J@2) 
DO 1112 1=1,42 
WRITE(2,1113) (D(1. J), J=15N2) 


1112 CONTINUE 
1113 FORMAT (10F8. 2) 


RETURN 
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APPENDIX J: SAMPLE FILES--GRID INTERPOLATION EXAMPLE 


35 30) «(50 40 10.0 0.5 


Figure J1. INTPGRD 


#1D MODS 
#D PARAM.3,PARAM.5 
PARAMETER(10=35, JO=30) 
PARAMETER ( 1Q2=50 , JQ2=40) 
PARAMETER(KO=2000) 
#D INTERPO.18 
11 FORMAT(15F5.0) 


Figure J2. INTPUPD 
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Figure J3. 
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APPENDIX K: NOTATION 


a or a(x,y) 


ol S&S Pl 


ow 


e or c(x,y) 


Cy or cy (x,y) 


Wave amplitude function 

Function of the bottom slope 
Function of dependent variables 
Function of the bottom slope 
Function of dependent variables 
Wave celerity 

Group velocity 

Deepwater wave celerity 

Energy dissipation function 
Energy dissipation function 

Wave energy 

Function of dependent variables 
Function of dependent variables 
Gravitational acceleration constant 
Function of dependent variables 
Function of dependent variables 
Water depth 

Water depth at incipient breaking 
Wave height 

Wave height at incipient breaking 
Deepwater wave height 

Arbitrary subscript designating the x-direction 


Complex number equal to v-1 


Unit vector in the x-direction 
Specific subscript designating the x-direction 


Arbitrary subscript designating the y-direction 


Unit vector in the y-direction 

Specific subscript designating the y-direction 
Wave number 

Wave length at incipient breaking 

Deepwater wave length 

Bottom slope 

Total number of grid cells in the x-direction 
Total number of grid cells in the y-direction 


Flow across a hydraulic jump 


K2 


s or s(x,y) Wave phase function 


s Subscript denoting stable wave conditions 
T Wave period 
W Weighting factor 
x Coordinate direction 
x! Coordinate direction 
y Coordinate direction 
ir Coordinate direction 
Yy Water depth on the high end of a hydraulic jump 
Yo Water depth on the low end of a hydraulic jump 
a Weighting factor 
a* Coefficient 
8 Coefficient 
Y Coefficient 


Rate of energy loss 


AX Grid size in the x-direction 
Ay Grid size in the y-direction 
i) Wave angle 
8, Contour angle 

8, Deepwater wave angle 

K Rate of energy dissipation coefficient 
ka Coefficient 
Ke Refraction coefficient 
Ky Shoaling coefficient 
Ws Coefficient 

p Water density 

(oj Angular wave frequency 

C0) Velocity potential function 


Mathematical symbols 


f) Partial differentiation 

V Horizontal gradient operator 
Vector dot product 

x Vector cross product 


| | Absolute value 


K3 


